Large Differences in Herbivore Performance Emerge from Simple Herbivore Behaviors and Fine-Scale Spatial Heterogeneity in Phytochemistry
Vincent S. Pan, Enakshi Ghosh, Paul J. Ode, William C. Wetzel, Kadeem J. Gilbert, Ian S. Pearse
July 11, 2024
1 Setup
1.1 Load libraries
library(survival)
library(piecewiseSEM)
source("spat1f/init_analysis.R")Data last updated on 2024-05-031.2 Handy functions not in spat1f
trend_by_pre_wt <- function(model, term){
emmeans::emtrends(model, var = "cat_pre_wt_log_scale", specs = term) %>%
pairs()
}
contrast_by_pre_wt <- function(model, term, type = c("pairwise", "trt.vs.ctrl"), at = c(-2, 2)){
type <- match.arg(type)
f <- as.formula(sprintf("%s ~ %s | cat_pre_wt_log_scale", type, term))
emmeans::emmeans(model,
f,
var = term,
at = list(cat_pre_wt_log_scale = at))
}
check_mod <- function(model){
DHARMa::simulateResiduals(model) %>%
plot()
performance::posterior_predictive_check(model)
}
# Get posterior draws
epred_draws <- function(model, newdata){
res <- brms::posterior_epred(model,
newdata = newdata,
allow_new_levels = TRUE,
re_formula = NA) %>%
t()
colnames(res) <- paste0("draw_", seq_len(ncol(res)))
cbind(newdata,res) %>%
gather(key = draw, value = val, contains("draw_"))
}1.3 Quick cleaning
# Main data.frame for analyses
d <- ref_data %>%
filter(error == 0) %>%
mutate(
cat_pre_wt_log = log(cat_pre_wt),
cat_pre_wt_log_scale = as.numeric(scale(log(cat_pre_wt))),
cat_size = ifelse(cat_pre_wt < median(cat_pre_wt), "small", "big"),
has_var = ifelse(var_trt != "constant", 1, 0),
mean_toxic_conc_scale = as.numeric(scale(mean_toxic_conc)),
mean_toxic_scale = as.numeric(scale(mean_toxic)),
area_herb_log_scale = as.numeric(scale(log(area_herb+1))),
var_toxic_12_scale = as.numeric(scale(var_toxic_12)),
on_toxic_logit_scale = as.numeric(scale(adjust_prop(on_toxic,
trans = "emp",
nudge.method = "none",
na.action = "ignore"))),
sl_mean_obs_log_scale = as.numeric(scale(log(sl_mean_obs))),
cat_pre_wt_log_scale_sq = as.numeric(scale(log(cat_pre_wt)^2)),
RGR_scale = as.numeric(scale(RGR)),
sl_mean_pred1 = shape_1 * scale_1, # Exploration state & on more toxic diet
sl_mean_pred2 = shape_2 * scale_2, # Resting/feeding state & on more toxic diet
sl_mean_pred3 = shape_3 * scale_3, # Exploration state & on less toxic diet
sl_mean_pred4 = shape_4 * scale_4, # Resting/feeding state & on less toxic diet
k1 = scale_1 / scale_3, # state 1 (exploration)
k2 = scale_2 / scale_4, # state 2 (resting/feeding)
prop_explore_logit = qlogis(prop_explore),
prop_explore_logit_scale = as.numeric(scale(prop_explore_logit))
)
# Subset of data.frame for SEM
d2 <- d %>%
filter(
var_trt != "constant"
) %>%
filter(
!is.na(area_herb_log_scale) &
!is.na(var_toxic_12_scale) &
!is.na(mean_toxic_conc_scale) &
!is.na(sl_mean_obs) &
!is.na(prop_explore) &
!is.na(RGR)) %>%
mutate(
var_high = ifelse(var_trt == "high_var", 1, 0),
beta_red = ifelse(as.character(beta) == 5, 1, 0),
beta_white = ifelse(as.character(beta) == 0, 1, 0),
beta_red_cat_pre_wt_log_scale = beta_red * cat_pre_wt_log_scale,
beta_white_cat_pre_wt_log_scale = beta_white * cat_pre_wt_log_scale,
var_high_cat_pre_wt_log_scale = var_high * cat_pre_wt_log_scale,
beta_numeric_scale = as.numeric(scale(as.numeric(beta)))
)
# Subset of data.frame for more detailed analyses
d3 <- d %>%
filter(
var_trt != "constant"
) %>%
filter_at(vars(contains("shape_"), contains("scale_[0-9]"), contains("sl_mean_pred")),
all_vars(. > 0)) %>%
mutate(
var_high = ifelse(var_trt == "high_var", 1, 0),
beta_red = ifelse(as.character(beta) == 5, 1, 0),
beta_white = ifelse(as.character(beta) == 0, 1, 0),
beta_red_cat_pre_wt_log_scale = beta_red * cat_pre_wt_log_scale,
beta_white_cat_pre_wt_log_scale = beta_white * cat_pre_wt_log_scale,
var_high_cat_pre_wt_log_scale = var_high * cat_pre_wt_log_scale,
beta_numeric_scale = as.numeric(scale(as.numeric(beta)))
)2 Herbivore Performance
2.1 RGR
2.1.1 Model summary
rgr_m <- glmmTMB(
RGR ~
(var_trt + beta) * cat_pre_wt_log_scale +
I(cat_pre_wt_log_scale^2) + # residual plot show significant non-linearity
(1|session_id),
data = d %>%
filter(
var_trt != "constant"
)
); summary(rgr_m) Family: gaussian ( identity )
Formula:
RGR ~ (var_trt + beta) * cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) +
(1 | session_id)
Data: d %>% filter(var_trt != "constant")
AIC BIC logLik deviance df.resid
-879.2 -849.4 450.6 -901.2 100
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 3.979e-06 0.001995
Residual 1.603e-05 0.004004
Number of obs: 111, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 1.6e-05
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.0162504 0.0013457 12.075 < 2e-16 ***
var_trtlow_var -0.0034872 0.0007740 -4.506 6.62e-06 ***
beta0 -0.0003065 0.0009664 -0.317 0.75112
beta5 -0.0013545 0.0009173 -1.477 0.13978
cat_pre_wt_log_scale -0.0048759 0.0011991 -4.066 4.78e-05 ***
I(cat_pre_wt_log_scale^2) -0.0021953 0.0006731 -3.262 0.00111 **
var_trtlow_var:cat_pre_wt_log_scale 0.0025873 0.0007897 3.276 0.00105 **
beta0:cat_pre_wt_log_scale 0.0011924 0.0009491 1.256 0.20897
beta5:cat_pre_wt_log_scale 0.0024731 0.0009710 2.547 0.01087 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(rgr_m, type = "III")Analysis of Deviance Table (Type III Wald chisquare tests)
Response: RGR
Chisq Df Pr(>Chisq)
(Intercept) 145.8141 1 < 2.2e-16 ***
var_trt 20.3008 1 6.617e-06 ***
beta 2.3963 2 0.301758
cat_pre_wt_log_scale 16.5348 1 4.777e-05 ***
I(cat_pre_wt_log_scale^2) 10.6374 1 0.001108 **
var_trt:cat_pre_wt_log_scale 10.7334 1 0.001052 **
beta:cat_pre_wt_log_scale 6.4936 2 0.038899 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 12.1.2 Post-hoc contrasts
trend_by_pre_wt(rgr_m, "var_trt") contrast estimate SE df t.ratio p.value
high_var - low_var -0.00259 0.00079 100 -3.276 0.0014
Results are averaged over the levels of: beta contrast_by_pre_wt(rgr_m, "var_trt")$emmeans
cat_pre_wt_log_scale = -2:
var_trt emmean SE df lower.CL upper.CL
high_var 0.014224 0.00355 100 0.007178 0.02127
low_var 0.005562 0.00312 100 -0.000635 0.01176
cat_pre_wt_log_scale = 2:
var_trt emmean SE df lower.CL upper.CL
high_var -0.000393 0.00208 100 -0.004512 0.00373
low_var 0.001295 0.00238 100 -0.003421 0.00601
Results are averaged over the levels of: beta
Confidence level used: 0.95
$contrasts
cat_pre_wt_log_scale = -2:
contrast estimate SE df t.ratio p.value
high_var - low_var 0.00866 0.00181 100 4.790 <.0001
cat_pre_wt_log_scale = 2:
contrast estimate SE df t.ratio p.value
high_var - low_var -0.00169 0.00171 100 -0.988 0.3256
Results are averaged over the levels of: beta trend_by_pre_wt(rgr_m, "beta") contrast estimate SE df t.ratio p.value
(beta-5) - beta0 -0.00119 0.000949 100 -1.256 0.4232
(beta-5) - beta5 -0.00247 0.000971 100 -2.547 0.0329
beta0 - beta5 -0.00128 0.000940 100 -1.363 0.3644
Results are averaged over the levels of: var_trt
P value adjustment: tukey method for comparing a family of 3 estimates contrast_by_pre_wt(rgr_m, "beta")$emmeans
cat_pre_wt_log_scale = -2:
beta emmean SE df lower.CL upper.CL
-5 0.012890 0.00382 100 0.005313 0.02047
0 0.010199 0.00333 100 0.003591 0.01681
5 0.006589 0.00319 100 0.000255 0.01292
cat_pre_wt_log_scale = 2:
beta emmean SE df lower.CL upper.CL
-5 -0.001439 0.00226 100 -0.005922 0.00304
0 0.000639 0.00231 100 -0.003940 0.00522
5 0.002153 0.00257 100 -0.002946 0.00725
Results are averaged over the levels of: var_trt
Confidence level used: 0.95
$contrasts
cat_pre_wt_log_scale = -2:
contrast estimate SE df t.ratio p.value
(beta-5) - beta0 0.00269 0.00221 100 1.218 0.4456
(beta-5) - beta5 0.00630 0.00223 100 2.820 0.0158
beta0 - beta5 0.00361 0.00211 100 1.711 0.2062
cat_pre_wt_log_scale = 2:
contrast estimate SE df t.ratio p.value
(beta-5) - beta0 -0.00208 0.00205 100 -1.015 0.5689
(beta-5) - beta5 -0.00359 0.00206 100 -1.746 0.1935
beta0 - beta5 -0.00151 0.00211 100 -0.719 0.7529
Results are averaged over the levels of: var_trt
P value adjustment: tukey method for comparing a family of 3 estimates 2.1.3 Misc stats
r2(rgr_m, tolerance = 10^-10)# R2 for Mixed Models
Conditional R2: 0.555
Marginal R2: 0.4442.1.4 Check model
check_mod(rgr_m)
2.1.5 Plot
g1 <- marginal_effects(rgr_m, terms = c("cat_pre_wt_log_scale","beta")) %>%
ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = yhat)) +
geom_ribbon(aes(ymax = upper, ymin = lower, fill = beta), alpha = 0.2) +
geom_line(aes(color = beta), linewidth = 2) +
geom_point(
data = rgr_m$frame,
aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = RGR, color = beta),
size = 3, alpha = 0.8
) +
theme_bw(base_size = 15) +
theme(legend.position = "top") +
scale_x_continuous(trans = "log10", labels = fancy_scientific) +
scale_color_brewer(type = "qual", aesthetics = c("fill","color")) +
labs(x = "Cat pre-weight (g)", y = expression(RGR~(hour^-1)),
color = expression(beta), fill = expression(beta))
g2 <- marginal_effects(rgr_m, terms = c("cat_pre_wt_log_scale","var_trt")) %>%
ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = yhat)) +
geom_ribbon(aes(ymax = upper, ymin = lower, fill = var_trt), alpha = 0.2) +
geom_line(aes(color = var_trt), linewidth = 2) +
geom_point(
data = rgr_m$frame,
aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = RGR, color = var_trt),
size = 3, alpha = 0.8
) +
theme_bw(base_size = 15) +
scale_x_continuous(trans = "log10", labels = fancy_scientific) +
theme(legend.position = "top") +
scale_color_discrete(type = c("navy","skyblue"),
aesthetics = c("fill","color"),
label = c("High", "Low")) +
labs(x = "Cat pre-weight (g)", y = expression(RGR~(hour^-1)),
color = "Variation", fill = "Variation")
ggarrange(g1, g2 + labs(y = ""))
2.2 Time to pupation
2.2.1 Model summary
pupt_m_full <- glmmTMB(
time_to_pupation ~
(var_trt + beta) * cat_pre_wt_log_scale +
I(cat_pre_wt_log_scale^2) + # Residuals show significant non-linearity
(1|session_id),
data = d %>%
filter(var_trt != "constant"),
family = Gamma(link = "log"),
); summary(pupt_m_full) Family: Gamma ( log )
Formula: time_to_pupation ~ (var_trt + beta) * cat_pre_wt_log_scale +
I(cat_pre_wt_log_scale^2) + (1 | session_id)
Data: d %>% filter(var_trt != "constant")
AIC BIC logLik deviance df.resid
338.1 366.3 -158.0 316.1 85
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 2.839e-11 5.328e-06
Number of obs: 96, groups: session_id, 5
Dispersion estimate for Gamma family (sigma^2): 0.0254
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.180302 0.036853 59.16 < 2e-16 ***
var_trtlow_var 0.102145 0.033998 3.00 0.00266 **
beta0 -0.001845 0.040823 -0.05 0.96395
beta5 0.059543 0.041181 1.45 0.14821
cat_pre_wt_log_scale -0.428662 0.032643 -13.13 < 2e-16 ***
I(cat_pre_wt_log_scale^2) -0.081267 0.019069 -4.26 2.03e-05 ***
var_trtlow_var:cat_pre_wt_log_scale 0.026518 0.038876 0.68 0.49516
beta0:cat_pre_wt_log_scale -0.030363 0.043195 -0.70 0.48210
beta5:cat_pre_wt_log_scale -0.126801 0.047095 -2.69 0.00709 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(pupt_m_full, type = "III")Analysis of Deviance Table (Type III Wald chisquare tests)
Response: time_to_pupation
Chisq Df Pr(>Chisq)
(Intercept) 3500.2207 1 < 2.2e-16 ***
var_trt 9.0264 1 0.002661 **
beta 2.7351 2 0.254736
cat_pre_wt_log_scale 172.4448 1 < 2.2e-16 ***
I(cat_pre_wt_log_scale^2) 18.1614 1 2.029e-05 ***
var_trt:cat_pre_wt_log_scale 0.4653 1 0.495162
beta:cat_pre_wt_log_scale 7.5589 2 0.022836 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1pupt_m <- glmmTMB(
time_to_pupation ~
var_trt + beta * cat_pre_wt_log_scale +
I(cat_pre_wt_log_scale^2) + # Residuals show significant non-linearity
(1|session_id),
data = d %>%
filter(var_trt != "constant"),
family = Gamma(link = "log"),
); summary(pupt_m) Family: Gamma ( log )
Formula:
time_to_pupation ~ var_trt + beta * cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) +
(1 | session_id)
Data: d %>% filter(var_trt != "constant")
AIC BIC logLik deviance df.resid
336.5 362.2 -158.3 316.5 86
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 2.159e-11 4.647e-06
Number of obs: 96, groups: session_id, 5
Dispersion estimate for Gamma family (sigma^2): 0.0255
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.180555 0.036937 59.03 < 2e-16 ***
var_trtlow_var 0.107447 0.033197 3.24 0.00121 **
beta0 -0.003578 0.040842 -0.09 0.93019
beta5 0.059536 0.041271 1.44 0.14915
cat_pre_wt_log_scale -0.420147 0.030221 -13.90 < 2e-16 ***
I(cat_pre_wt_log_scale^2) -0.083365 0.018881 -4.42 1.01e-05 ***
beta0:cat_pre_wt_log_scale -0.026011 0.042849 -0.61 0.54383
beta5:cat_pre_wt_log_scale -0.126834 0.047183 -2.69 0.00719 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(pupt_m, type = "III")Analysis of Deviance Table (Type III Wald chisquare tests)
Response: time_to_pupation
Chisq Df Pr(>Chisq)
(Intercept) 3485.0791 1 < 2.2e-16 ***
var_trt 10.4759 1 0.001209 **
beta 2.8053 2 0.245947
cat_pre_wt_log_scale 193.2743 1 < 2.2e-16 ***
I(cat_pre_wt_log_scale^2) 19.4946 1 1.009e-05 ***
beta:cat_pre_wt_log_scale 7.6890 2 0.021397 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 12.2.2 Post-hoc contrasts
contrast_by_pre_wt(pupt_m, "var_trt")$emmeans
cat_pre_wt_log_scale = -2:
var_trt emmean SE df asymp.LCL asymp.UCL
high_var 2.808 0.0769 Inf 2.657 2.96
low_var 2.915 0.0818 Inf 2.755 3.08
cat_pre_wt_log_scale = 2:
var_trt emmean SE df asymp.LCL asymp.UCL
high_var 0.924 0.0672 Inf 0.792 1.06
low_var 1.031 0.0725 Inf 0.889 1.17
Results are averaged over the levels of: beta
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
cat_pre_wt_log_scale = -2:
contrast estimate SE df z.ratio p.value
high_var - low_var -0.107 0.0332 Inf -3.237 0.0012
cat_pre_wt_log_scale = 2:
contrast estimate SE df z.ratio p.value
high_var - low_var -0.107 0.0332 Inf -3.237 0.0012
Results are averaged over the levels of: beta
Results are given on the log (not the response) scale. trend_by_pre_wt(pupt_m, "beta") contrast estimate SE df z.ratio p.value
(beta-5) - beta0 0.026 0.0428 Inf 0.607 0.8163
(beta-5) - beta5 0.127 0.0472 Inf 2.688 0.0197
beta0 - beta5 0.101 0.0473 Inf 2.132 0.0834
Results are averaged over the levels of: var_trt
P value adjustment: tukey method for comparing a family of 3 estimates contrast_by_pre_wt(pupt_m, "beta")$emmeans
cat_pre_wt_log_scale = -2:
beta emmean SE df asymp.LCL asymp.UCL
-5 2.741 0.0941 Inf 2.557 2.93
0 2.790 0.0971 Inf 2.599 2.98
5 3.054 0.1088 Inf 2.841 3.27
cat_pre_wt_log_scale = 2:
beta emmean SE df asymp.LCL asymp.UCL
-5 1.061 0.0845 Inf 0.895 1.23
0 1.005 0.0805 Inf 0.847 1.16
5 0.866 0.0917 Inf 0.687 1.05
Results are averaged over the levels of: var_trt
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
cat_pre_wt_log_scale = -2:
contrast estimate SE df z.ratio p.value
(beta-5) - beta0 -0.0484 0.1018 Inf -0.476 0.8828
(beta-5) - beta5 -0.3132 0.1133 Inf -2.763 0.0158
beta0 - beta5 -0.2648 0.1136 Inf -2.331 0.0516
cat_pre_wt_log_scale = 2:
contrast estimate SE df z.ratio p.value
(beta-5) - beta0 0.0556 0.0875 Inf 0.636 0.8005
(beta-5) - beta5 0.1941 0.0915 Inf 2.122 0.0854
beta0 - beta5 0.1385 0.0928 Inf 1.492 0.2945
Results are averaged over the levels of: var_trt
Results are given on the log (not the response) scale.
P value adjustment: tukey method for comparing a family of 3 estimates 2.2.3 Misc stats
r2(pupt_m, tolerance = 10^-10)Random effect variances not available. Returned R2 does not account for random effects.# R2 for Mixed Models
Conditional R2: NA
Marginal R2: 0.8812.2.4 Check model
check_mod(pupt_m)
2.2.5 Plot
g3 <- marginal_effects(pupt_m, terms = c("cat_pre_wt_log_scale","beta")) %>%
ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = yhat)) +
geom_ribbon(aes(ymax = upper, ymin = lower, fill = beta), alpha = 0.2) +
geom_line(aes(color = beta), linewidth = 2) +
geom_point(
data = pupt_m$frame,
aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = time_to_pupation, color = beta),
size = 3, alpha = 0.8
) +
theme_bw(base_size = 15) +
theme(legend.position = "top") +
scale_y_continuous(trans = "log2") +
scale_x_continuous(trans = "log10", labels = fancy_scientific) +
scale_color_brewer(type = "qual", aesthetics = c("fill","color")) +
labs(x = "Cat pre-weight (g)", y = "Time to pupation (days)",
color = expression(beta), fill = expression(beta))
g4 <- marginal_effects(pupt_m, terms = c("var_trt")) %>%
ggplot(aes(x = var_trt, y = yhat)) +
geom_point(
data = pupt_m$frame,
aes(x = var_trt, y = time_to_pupation, color = var_trt),
position = position_jitter(width = 0.2),
size = 3, alpha = 0.8
) +
geom_pointrange(
aes(ymax = upper, ymin = lower),
color = "black",
size = 1, linewidth = 2, shape = 1,
) +
theme_bw(base_size = 15) +
scale_y_continuous(trans = "log2") +
theme(legend.position = "top") +
scale_color_discrete(type = c("navy","skyblue"),
aesthetics = c("color"),
label = c("High", "Low")) +
labs(x = "Variation", y = "Time to pupation (days)",
color = "Variation")
ggarrange(g3, g4 + labs(y = ""))
2.3 Eclosure
2.3.1 Model summary
eclosure_m_full <- glmmTMB(
eclosed ~
(var_trt + beta) * cat_pre_wt_log_scale +
(1|session_id),
data = d %>%
filter(var_trt != "constant"),
family = binomial(link = "logit"),
); summary(eclosure_m_full) Family: binomial ( logit )
Formula:
eclosed ~ (var_trt + beta) * cat_pre_wt_log_scale + (1 | session_id)
Data: d %>% filter(var_trt != "constant")
AIC BIC logLik deviance df.resid
185.3 210.9 -83.6 167.3 119
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 0.1661 0.4076
Number of obs: 128, groups: session_id, 5
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.53530 0.41050 1.304 0.192
var_trtlow_var 0.01707 0.37507 0.046 0.964
beta0 -0.21544 0.46265 -0.466 0.641
beta5 -0.50749 0.45596 -1.113 0.266
cat_pre_wt_log_scale 0.10452 0.36089 0.290 0.772
var_trtlow_var:cat_pre_wt_log_scale 0.17526 0.36489 0.480 0.631
beta0:cat_pre_wt_log_scale 0.24778 0.43553 0.569 0.569
beta5:cat_pre_wt_log_scale 0.20159 0.44699 0.451 0.652car::Anova(eclosure_m_full, type = "III")Analysis of Deviance Table (Type III Wald chisquare tests)
Response: eclosed
Chisq Df Pr(>Chisq)
(Intercept) 1.7004 1 0.1922
var_trt 0.0021 1 0.9637
beta 1.2502 2 0.5352
cat_pre_wt_log_scale 0.0839 1 0.7721
var_trt:cat_pre_wt_log_scale 0.2307 1 0.6310
beta:cat_pre_wt_log_scale 0.3659 2 0.8328eclosure_m <- glmmTMB(
eclosed ~
(var_trt + beta) + cat_pre_wt_log_scale +
(1|session_id),
data = d %>%
filter(var_trt != "constant"),
family = binomial(link = "logit"),
); summary(eclosure_m) Family: binomial ( logit )
Formula:
eclosed ~ (var_trt + beta) + cat_pre_wt_log_scale + (1 | session_id)
Data: d %>% filter(var_trt != "constant")
AIC BIC logLik deviance df.resid
179.9 197.0 -84.0 167.9 122
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 0.1676 0.4094
Number of obs: 128, groups: session_id, 5
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.551286 0.413578 1.333 0.183
var_trtlow_var 0.005064 0.373402 0.014 0.989
beta0 -0.237692 0.460689 -0.516 0.606
beta5 -0.523521 0.456587 -1.147 0.252
cat_pre_wt_log_scale 0.330169 0.220049 1.500 0.134car::Anova(eclosure_m, type = "II")Analysis of Deviance Table (Type II Wald chisquare tests)
Response: eclosed
Chisq Df Pr(>Chisq)
var_trt 0.0002 1 0.9892
beta 1.3212 2 0.5165
cat_pre_wt_log_scale 2.2513 1 0.13352.3.2 Post-hoc contrasts
2.3.3 Misc stats
r2(eclosure_m, tolerance = 10^-10)# R2 for Mixed Models
Conditional R2: 0.092
Marginal R2: 0.0462.3.4 Check model
check_mod(eclosure_m)
2.3.5 Plot
2.4 Survival
2.4.1 Model summary
surv_m_full <- coxph(
Surv(surv_time, observed_dead) ~
(var_trt + beta) * cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) +
strata(session_id),
data = d %>%
filter(var_trt != "constant") %>%
mutate(
beta = as.factor(as.character(beta))
),
); summary(surv_m_full)Call:
coxph(formula = Surv(surv_time, observed_dead) ~ (var_trt + beta) *
cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) + strata(session_id),
data = d %>% filter(var_trt != "constant") %>% mutate(beta = as.factor(as.character(beta))))
n= 128, number of events= 105
coef exp(coef) se(coef) z Pr(>|z|)
var_trtlow_var -0.1115 0.8945 0.2124 -0.525 0.599727
beta0 -0.0896 0.9143 0.2579 -0.347 0.728306
beta5 -0.2521 0.7772 0.2605 -0.968 0.333144
cat_pre_wt_log_scale 0.9641 2.6224 0.2556 3.772 0.000162
I(cat_pre_wt_log_scale^2) 0.2450 1.2776 0.1447 1.693 0.090408
var_trtlow_var:cat_pre_wt_log_scale -0.1626 0.8499 0.1998 -0.814 0.415661
beta0:cat_pre_wt_log_scale -0.2230 0.8001 0.2309 -0.966 0.334053
beta5:cat_pre_wt_log_scale -0.4789 0.6195 0.2584 -1.854 0.063807
var_trtlow_var
beta0
beta5
cat_pre_wt_log_scale ***
I(cat_pre_wt_log_scale^2) .
var_trtlow_var:cat_pre_wt_log_scale
beta0:cat_pre_wt_log_scale
beta5:cat_pre_wt_log_scale .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
var_trtlow_var 0.8945 1.1179 0.5899 1.356
beta0 0.9143 1.0937 0.5515 1.516
beta5 0.7772 1.2867 0.4665 1.295
cat_pre_wt_log_scale 2.6224 0.3813 1.5892 4.328
I(cat_pre_wt_log_scale^2) 1.2776 0.7827 0.9621 1.696
var_trtlow_var:cat_pre_wt_log_scale 0.8499 1.1766 0.5746 1.257
beta0:cat_pre_wt_log_scale 0.8001 1.2499 0.5089 1.258
beta5:cat_pre_wt_log_scale 0.6195 1.6142 0.3734 1.028
Concordance= 0.637 (se = 0.047 )
Likelihood ratio test= 22 on 8 df, p=0.005
Wald test = 22.31 on 8 df, p=0.004
Score (logrank) test = 24.57 on 8 df, p=0.002drop1(surv_m_full, test = "Chisq")Single term deletions
Model:
Surv(surv_time, observed_dead) ~ (var_trt + beta) * cat_pre_wt_log_scale +
I(cat_pre_wt_log_scale^2) + strata(session_id)
Df AIC LRT Pr(>Chi)
<none> 459.50
I(cat_pre_wt_log_scale^2) 1 460.40 2.90 0.08841 .
strata(session_id) 0 787.87 328.37
var_trt:cat_pre_wt_log_scale 1 458.16 0.67 0.41395
beta:cat_pre_wt_log_scale 2 458.98 3.48 0.17532
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1surv_m <- coxph(
Surv(surv_time, observed_dead) ~
(var_trt + beta) + cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) +
strata(session_id),
data = d %>%
filter(var_trt != "constant") %>%
mutate(
beta = as.factor(as.character(beta))
),
); summary(surv_m)Call:
coxph(formula = Surv(surv_time, observed_dead) ~ (var_trt + beta) +
cat_pre_wt_log_scale + I(cat_pre_wt_log_scale^2) + strata(session_id),
data = d %>% filter(var_trt != "constant") %>% mutate(beta = as.factor(as.character(beta))))
n= 128, number of events= 105
coef exp(coef) se(coef) z Pr(>|z|)
var_trtlow_var -0.07077 0.93167 0.20845 -0.340 0.734212
beta0 -0.07410 0.92858 0.25808 -0.287 0.774026
beta5 -0.17722 0.83759 0.25804 -0.687 0.492206
cat_pre_wt_log_scale 0.65947 1.93376 0.18164 3.631 0.000283 ***
I(cat_pre_wt_log_scale^2) 0.29584 1.34425 0.14124 2.095 0.036205 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
var_trtlow_var 0.9317 1.0733 0.6192 1.402
beta0 0.9286 1.0769 0.5599 1.540
beta5 0.8376 1.1939 0.5051 1.389
cat_pre_wt_log_scale 1.9338 0.5171 1.3545 2.761
I(cat_pre_wt_log_scale^2) 1.3443 0.7439 1.0192 1.773
Concordance= 0.602 (se = 0.047 )
Likelihood ratio test= 18.08 on 5 df, p=0.003
Wald test = 18.29 on 5 df, p=0.003
Score (logrank) test = 19.42 on 5 df, p=0.002drop1(surv_m, test = "Chisq")Single term deletions
Model:
Surv(surv_time, observed_dead) ~ (var_trt + beta) + cat_pre_wt_log_scale +
I(cat_pre_wt_log_scale^2) + strata(session_id)
Df AIC LRT Pr(>Chi)
<none> 457.42
var_trt 1 455.53 0.12 0.7341150
beta 2 453.89 0.48 0.7883535
cat_pre_wt_log_scale 1 469.05 13.63 0.0002226 ***
I(cat_pre_wt_log_scale^2) 1 459.84 4.42 0.0355376 *
strata(session_id) 0 784.70 327.28
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 12.4.2 Post-hoc contrasts
2.4.3 Misc stats
r2(surv_m, tolerance = 10^-10) Nagelkerke's R2: 0.1352.4.4 Check model
cox.zph(surv_m) chisq df p
var_trt 1.60 1 0.21
beta 1.43 2 0.49
cat_pre_wt_log_scale 2.61 1 0.11
I(cat_pre_wt_log_scale^2) 1.05 1 0.31
GLOBAL 5.53 5 0.352.4.5 Plot
3 Structural Equation Model
3.1 Sub-models
3.1.1 RGR
subm_rgr <- glmmTMB(
RGR_scale ~
cat_pre_wt_log_scale +
cat_pre_wt_log_scale_sq +
sl_mean_obs_log_scale +
prop_explore_logit_scale *
var_toxic_12_scale +
mean_toxic_conc_scale +
area_herb_log_scale +
beta_numeric_scale +
(1|session_id),
data = d2
); summary(subm_rgr) Family: gaussian ( identity )
Formula: RGR_scale ~ cat_pre_wt_log_scale + cat_pre_wt_log_scale_sq +
sl_mean_obs_log_scale + prop_explore_logit_scale * var_toxic_12_scale +
mean_toxic_conc_scale + area_herb_log_scale + beta_numeric_scale +
(1 | session_id)
Data: d2
AIC BIC logLik deviance df.resid
119.5 149.0 -47.7 95.5 75
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 0.01245 0.1116
Residual 0.16744 0.4092
Number of obs: 87, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.167
Conditional model:
Estimate Std. Error z value
(Intercept) -0.02766 0.07257 -0.381
cat_pre_wt_log_scale -3.64523 0.41435 -8.797
cat_pre_wt_log_scale_sq -3.01946 0.45628 -6.618
sl_mean_obs_log_scale -0.13590 0.05626 -2.416
prop_explore_logit_scale 0.02370 0.04889 0.485
var_toxic_12_scale -0.01949 0.05115 -0.381
mean_toxic_conc_scale -0.25503 0.07147 -3.568
area_herb_log_scale 0.76522 0.09726 7.868
beta_numeric_scale -0.11265 0.04837 -2.329
prop_explore_logit_scale:var_toxic_12_scale 0.18687 0.04488 4.164
Pr(>|z|)
(Intercept) 0.703071
cat_pre_wt_log_scale < 2e-16 ***
cat_pre_wt_log_scale_sq 3.65e-11 ***
sl_mean_obs_log_scale 0.015710 *
prop_explore_logit_scale 0.627788
var_toxic_12_scale 0.703095
mean_toxic_conc_scale 0.000359 ***
area_herb_log_scale 3.61e-15 ***
beta_numeric_scale 0.019873 *
prop_explore_logit_scale:var_toxic_12_scale 3.13e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.1.2 Variance in toxin (12 hour)
subm_var_toxic <- glmmTMB(
var_toxic_12_scale ~
sl_mean_obs_log_scale +
beta_numeric_scale +
var_high +
(1|session_id),
data = d2
); summary(subm_var_toxic) Family: gaussian ( identity )
Formula:
var_toxic_12_scale ~ sl_mean_obs_log_scale + beta_numeric_scale +
var_high + (1 | session_id)
Data: d2
AIC BIC logLik deviance df.resid
178.2 193.0 -83.1 166.2 81
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 6.229e-11 7.892e-06
Residual 3.956e-01 6.290e-01
Number of obs: 87, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.396
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.75075 0.09477 -7.922 2.34e-15 ***
sl_mean_obs_log_scale 0.21267 0.07751 2.744 0.00607 **
beta_numeric_scale -0.15689 0.06811 -2.303 0.02126 *
var_high 1.56021 0.13868 11.250 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.1.3 Mean step length
subm_sl <- glmmTMB(
sl_mean_obs_log_scale ~
prop_explore_logit_scale +
cat_pre_wt_log_scale +
cat_pre_wt_log_scale_sq +
on_toxic_logit_scale +
(1|session_id),
data = d2,
); summary(subm_sl) Family: gaussian ( identity )
Formula:
sl_mean_obs_log_scale ~ prop_explore_logit_scale + cat_pre_wt_log_scale +
cat_pre_wt_log_scale_sq + on_toxic_logit_scale + (1 | session_id)
Data: d2
AIC BIC logLik deviance df.resid
223.4 240.7 -104.7 209.4 80
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 4.610e-10 2.147e-05
Residual 6.502e-01 8.064e-01
Number of obs: 87, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.65
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.06163 0.09646 0.639 0.52287
prop_explore_logit_scale -0.26175 0.08987 -2.912 0.00359 **
cat_pre_wt_log_scale 1.37948 0.70121 1.967 0.04915 *
cat_pre_wt_log_scale_sq 1.61095 0.75569 2.132 0.03303 *
on_toxic_logit_scale 0.21435 0.09489 2.259 0.02389 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.1.4 Mean toxin ingested
subm_toxin_ingested <- glmmTMB(
mean_toxic_conc_scale ~
beta_numeric_scale +
on_toxic_logit_scale +
var_high +
(1|session_id),
data = d2
); summary(subm_toxin_ingested) Family: gaussian ( identity )
Formula:
mean_toxic_conc_scale ~ beta_numeric_scale + on_toxic_logit_scale +
var_high + (1 | session_id)
Data: d2
AIC BIC logLik deviance df.resid
145.7 160.5 -66.8 133.7 81
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 6.461e-10 2.542e-05
Residual 2.722e-01 5.217e-01
Number of obs: 87, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.272
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.18349 0.07861 2.334 0.0196 *
beta_numeric_scale -0.11683 0.05808 -2.012 0.0443 *
on_toxic_logit_scale 0.28426 0.06393 4.447 8.72e-06 ***
var_high -0.80646 0.11613 -6.945 3.79e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.1.5 Mean time on toxin
subm_on_toxic <- glmmTMB(
on_toxic_logit_scale ~
var_high + beta_numeric_scale * cat_pre_wt_log_scale +
(1|session_id),
family = gaussian(),
data = d2
); summary(subm_on_toxic) Family: gaussian ( identity )
Formula:
on_toxic_logit_scale ~ var_high + beta_numeric_scale * cat_pre_wt_log_scale +
(1 | session_id)
Data: d2
AIC BIC logLik deviance df.resid
230.0 247.3 -108.0 216.0 80
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 0.01667 0.1291
Residual 0.68757 0.8292
Number of obs: 87, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.688
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.22509 0.14140 1.592 0.1114
var_high -0.41809 0.18032 -2.319 0.0204 *
beta_numeric_scale -0.12409 0.09766 -1.271 0.2039
cat_pre_wt_log_scale -0.19306 0.13787 -1.400 0.1614
beta_numeric_scale:cat_pre_wt_log_scale -0.23850 0.11732 -2.033 0.0421 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.1.6 Area consumption
subm_herb <- glmmTMB(
area_herb_log_scale ~
var_high * cat_pre_wt_log_scale +
(1|session_id),
data = d2,
); summary(subm_herb) Family: gaussian ( identity )
Formula: area_herb_log_scale ~ var_high * cat_pre_wt_log_scale + (1 |
session_id)
Data: d2
AIC BIC logLik deviance df.resid
117.9 132.7 -53.0 105.9 81
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 3.501e-11 5.917e-06
Residual 1.978e-01 4.448e-01
Number of obs: 87, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.198
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.04339 0.06963 0.623 0.53313
var_high 0.27072 0.10127 2.673 0.00751 **
cat_pre_wt_log_scale 0.55603 0.08378 6.637 3.2e-11 ***
var_high:cat_pre_wt_log_scale -0.26289 0.11558 -2.275 0.02293 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.1.7 Exploration
subm_exp <- glmmTMB(
prop_explore_logit_scale ~
var_high + beta_numeric_scale * cat_pre_wt_log_scale +
(1|session_id),
family = gaussian(),
data = d2
); summary(subm_exp) Family: gaussian ( identity )
Formula:
prop_explore_logit_scale ~ var_high + beta_numeric_scale * cat_pre_wt_log_scale +
(1 | session_id)
Data: d2
AIC BIC logLik deviance df.resid
250.1 267.3 -118.0 236.1 80
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 3.781e-10 1.944e-05
Residual 8.830e-01 9.397e-01
Number of obs: 87, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.883
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.17231 0.14397 -1.197 0.2314
var_high 0.15661 0.20355 0.769 0.4416
beta_numeric_scale -0.14059 0.11040 -1.273 0.2029
cat_pre_wt_log_scale -0.04358 0.12362 -0.353 0.7244
beta_numeric_scale:cat_pre_wt_log_scale 0.25977 0.13242 1.962 0.0498 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.2 d-separation test
sem_fit <- psem(
subm_rgr,
subm_herb,
subm_var_toxic,
subm_on_toxic,
subm_toxin_ingested,
subm_exp,
subm_sl,
var_toxic_12_scale %~~% on_toxic_logit_scale, # correlated errors
data = d2
)
sem_summary <- suppressWarnings(suppressMessages(
summary_psem(sem_fit,
no_standardize_x = c(),
.progressBar = FALSE,
direction = "var_toxic_12_scale -> area_herb_log_scale")
))Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.
Random effect variances not available. Returned R2 does not account for random effects.sem_summary$Cstat Fisher.C df P.Value
1 42.469 42 0.451knitr::kable(sem_summary$dTable)| Independ.Claim | Test.Type | DF | Crit.Value | P.Value | |
|---|---|---|---|---|---|
| mean_toxic_conc_scale ~ cat_pre_wt_log_scale + … | coef | 87 | -1.5712 | 0.1161 | |
| var_toxic_12_scale ~ cat_pre_wt_log_scale + … | coef | 87 | -0.1892 | 0.8499 | |
| prop_explore_logit_scale ~ cat_pre_wt_log_scale_sq + … | coef | 87 | 0.8624 | 0.3884 | |
| mean_toxic_conc_scale ~ cat_pre_wt_log_scale_sq + … | coef | 87 | 1.8002 | 0.0718 | |
| area_herb_log_scale ~ cat_pre_wt_log_scale_sq + … | coef | 87 | -0.7057 | 0.4804 | |
| on_toxic_logit_scale ~ cat_pre_wt_log_scale_sq + … | coef | 87 | 0.7944 | 0.4269 | |
| var_toxic_12_scale ~ cat_pre_wt_log_scale_sq + … | coef | 87 | 0.3081 | 0.7580 | |
| sl_mean_obs_log_scale ~ beta_numeric_scale + … | coef | 87 | 1.1020 | 0.2704 | |
| area_herb_log_scale ~ beta_numeric_scale + … | coef | 87 | 0.0200 | 0.9840 | |
| RGR_scale ~ var_high + … | coef | 87 | 0.6991 | 0.4845 | |
| sl_mean_obs_log_scale ~ var_high + … | coef | 87 | -1.5883 | 0.1122 | |
| on_toxic_logit_scale ~ RGR_scale + … | coef | 87 | -1.3664 | 0.1718 | |
| mean_toxic_conc_scale ~ sl_mean_obs_log_scale + … | coef | 87 | 0.7459 | 0.4557 | |
| area_herb_log_scale ~ sl_mean_obs_log_scale + … | coef | 87 | 0.8564 | 0.3918 | |
| mean_toxic_conc_scale ~ prop_explore_logit_scale + … | coef | 87 | 1.0591 | 0.2895 | |
| area_herb_log_scale ~ prop_explore_logit_scale + … | coef | 87 | 0.3920 | 0.6951 | |
| on_toxic_logit_scale ~ prop_explore_logit_scale + … | coef | 87 | 0.2361 | 0.8134 | |
| var_toxic_12_scale ~ prop_explore_logit_scale + … | coef | 87 | 1.1110 | 0.2666 | |
| area_herb_log_scale ~ mean_toxic_conc_scale + … | coef | 87 | -0.2457 | 0.8059 | |
| var_toxic_12_scale ~ mean_toxic_conc_scale + … | coef | 87 | -1.3086 | 0.1907 | |
| on_toxic_logit_scale ~ area_herb_log_scale + … | coef | 87 | -0.7213 | 0.4707 |
3.3 Model summary
3.3.1 Model stats
sem_summary$ChiSqNULLsem_summary$AIC AIC AICc K n
1 1264.813 1276.432 51 87sem_summary$R2 Response family link method Marginal Conditional
1 RGR_scale gaussian identity none 0.75 0.77
2 area_herb_log_scale gaussian identity none 0.42 NA
3 var_toxic_12_scale gaussian identity none 0.60 NA
4 on_toxic_logit_scale gaussian identity none 0.19 0.21
5 mean_toxic_conc_scale gaussian identity none 0.55 NA
6 prop_explore_logit_scale gaussian identity none 0.06 NA
7 sl_mean_obs_log_scale gaussian identity none 0.19 NA3.3.2 Standardized coefficients
knitr::kable(
cbind(sem_summary$coefficients,"arrow_size"=(abs(as.numeric(sem_summary$coefficients$Std.Estimate))*20))
)| Response | Predictor | Estimate | Std.Error | DF | Crit.Value | P.Value | Std.Estimate | star | arrow_size |
|---|---|---|---|---|---|---|---|---|---|
| RGR_scale | cat_pre_wt_log_scale | -3.6452 | 0.4144 | 87 | -8.7974 | 0.0000 | -3.5981 | *** | 71.962 |
| RGR_scale | cat_pre_wt_log_scale_sq | -3.0195 | 0.4563 | 87 | -6.6176 | 0.0000 | -2.7708 | *** | 55.416 |
| RGR_scale | sl_mean_obs_log_scale | -0.1359 | 0.0563 | 87 | -2.4156 | 0.0157 | -0.1448 | * | 2.896 |
| RGR_scale | prop_explore_logit_scale | 0.0237 | 0.0489 | 87 | 0.4848 | 0.6278 | 0.0274 | 0.548 | |
| RGR_scale | var_toxic_12_scale | -0.0195 | 0.0511 | 87 | -0.3811 | 0.7031 | -0.0230 | 0.460 | |
| RGR_scale | mean_toxic_conc_scale | -0.2550 | 0.0715 | 87 | -3.5683 | 0.0004 | -0.2353 | *** | 4.706 |
| RGR_scale | area_herb_log_scale | 0.7652 | 0.0973 | 87 | 7.8679 | 0.0000 | 0.5338 | *** | 10.676 |
| RGR_scale | beta_numeric_scale | -0.1126 | 0.0484 | 87 | -2.3287 | 0.0199 | -0.1337 | * | 2.674 |
| RGR_scale | prop_explore_logit_scale:var_toxic_12_scale | 0.1869 | 0.0449 | 87 | 4.1640 | 0.0000 | 0.2280 | *** | 4.560 |
| area_herb_log_scale | var_high | 0.2707 | 0.1013 | 87 | 2.6733 | 0.0075 | 0.2315 | ** | 4.630 |
| area_herb_log_scale | cat_pre_wt_log_scale | 0.5560 | 0.0838 | 87 | 6.6371 | 0.0000 | 0.7868 | *** | 15.736 |
| area_herb_log_scale | var_high:cat_pre_wt_log_scale | -0.2629 | 0.1156 | 87 | -2.2745 | 0.0229 | -0.2795 | * | 5.590 |
| var_toxic_12_scale | sl_mean_obs_log_scale | 0.2127 | 0.0775 | 87 | 2.7437 | 0.0061 | 0.1918 | ** | 3.836 |
| var_toxic_12_scale | beta_numeric_scale | -0.1569 | 0.0681 | 87 | -2.3033 | 0.0213 | -0.1575 | * | 3.150 |
| var_toxic_12_scale | var_high | 1.5602 | 0.1387 | 87 | 11.2503 | 0.0000 | 0.7875 | *** | 15.750 |
| on_toxic_logit_scale | var_high | -0.4181 | 0.1803 | 87 | -2.3186 | 0.0204 | -0.2235 | * | 4.470 |
| on_toxic_logit_scale | beta_numeric_scale | -0.1241 | 0.0977 | 87 | -1.2706 | 0.2039 | -0.1320 | 2.640 | |
| on_toxic_logit_scale | cat_pre_wt_log_scale | -0.1931 | 0.1379 | 87 | -1.4003 | 0.1614 | -0.1708 | 3.416 | |
| on_toxic_logit_scale | beta_numeric_scale:cat_pre_wt_log_scale | -0.2385 | 0.1173 | 87 | -2.0329 | 0.0421 | -0.2116 | * | 4.232 |
| mean_toxic_conc_scale | beta_numeric_scale | -0.1168 | 0.0581 | 87 | -2.0116 | 0.0443 | -0.1503 | * | 3.006 |
| mean_toxic_conc_scale | on_toxic_logit_scale | 0.2843 | 0.0639 | 87 | 4.4467 | 0.0000 | 0.3438 | *** | 6.876 |
| mean_toxic_conc_scale | var_high | -0.8065 | 0.1161 | 87 | -6.9447 | 0.0000 | -0.5214 | *** | 10.428 |
| prop_explore_logit_scale | var_high | 0.1566 | 0.2035 | 87 | 0.7694 | 0.4416 | 0.0808 | 1.616 | |
| prop_explore_logit_scale | beta_numeric_scale | -0.1406 | 0.1104 | 87 | -1.2734 | 0.2029 | -0.1443 | 2.886 | |
| prop_explore_logit_scale | cat_pre_wt_log_scale | -0.0436 | 0.1236 | 87 | -0.3526 | 0.7244 | -0.0372 | 0.744 | |
| prop_explore_logit_scale | beta_numeric_scale:cat_pre_wt_log_scale | 0.2598 | 0.1324 | 87 | 1.9617 | 0.0498 | 0.2224 | * | 4.448 |
| sl_mean_obs_log_scale | prop_explore_logit_scale | -0.2618 | 0.0899 | 87 | -2.9124 | 0.0036 | -0.2840 | ** | 5.680 |
| sl_mean_obs_log_scale | cat_pre_wt_log_scale | 1.3795 | 0.7012 | 87 | 1.9673 | 0.0491 | 1.2778 | * | 25.556 |
| sl_mean_obs_log_scale | cat_pre_wt_log_scale_sq | 1.6109 | 0.7557 | 87 | 2.1318 | 0.0330 | 1.3873 | * | 27.746 |
| sl_mean_obs_log_scale | on_toxic_logit_scale | 0.2144 | 0.0949 | 87 | 2.2589 | 0.0239 | 0.2244 | * | 4.488 |
| ~~var_toxic_12_scale | ~~on_toxic_logit_scale | 0.5323 | - | 87 | 5.7626 | 0.0000 | 0.5323 | *** | 10.646 |
3.3.3 Post hoc contrast
contrast_by_pre_wt(subm_herb, "var_high", "trt")$emmeans
cat_pre_wt_log_scale = -2:
var_high emmean SE df lower.CL upper.CL
0 -1.069 0.200 81 -1.467 -0.670
1 -0.272 0.198 81 -0.666 0.122
cat_pre_wt_log_scale = 2:
var_high emmean SE df lower.CL upper.CL
0 1.155 0.161 81 0.836 1.475
1 0.900 0.150 81 0.603 1.198
Confidence level used: 0.95
$contrasts
cat_pre_wt_log_scale = -2:
contrast estimate SE df t.ratio p.value
var_high1 - var_high0 0.797 0.281 81 2.830 0.0059
cat_pre_wt_log_scale = 2:
contrast estimate SE df t.ratio p.value
var_high1 - var_high0 -0.255 0.219 81 -1.162 0.2485sd(subm_herb$frame$var_high)/sd(subm_herb$frame$area_herb_log_scale)[1] 0.8549509emmeans::emtrends(subm_exp,
~ beta_numeric_scale | cat_pre_wt_log_scale,
var = "beta_numeric_scale",
at = list(cat_pre_wt_log_scale = c(-2,2)))cat_pre_wt_log_scale = -2:
beta_numeric_scale beta_numeric_scale.trend SE df lower.CL upper.CL
7.88e-17 -0.660 0.323 80 -1.30 -0.0173
cat_pre_wt_log_scale = 2:
beta_numeric_scale beta_numeric_scale.trend SE df lower.CL upper.CL
7.88e-17 0.379 0.246 80 -0.11 0.8677
Results are averaged over the levels of: var_high
Confidence level used: 0.95 sd(subm_exp$frame$beta_numeric_scale) * sd(subm_exp$frame$prop_explore_logit_scale)[1] 0.974435emmeans::emtrends(subm_on_toxic,
~ beta_numeric_scale | cat_pre_wt_log_scale,
var = "beta_numeric_scale",
at = list(cat_pre_wt_log_scale = c(-2,2)))cat_pre_wt_log_scale = -2:
beta_numeric_scale beta_numeric_scale.trend SE df lower.CL upper.CL
7.88e-17 0.353 0.286 80 -0.217 0.923
cat_pre_wt_log_scale = 2:
beta_numeric_scale beta_numeric_scale.trend SE df lower.CL upper.CL
7.88e-17 -0.601 0.217 80 -1.034 -0.169
Results are averaged over the levels of: var_high
Confidence level used: 0.95 sd(subm_on_toxic$frame$beta_numeric_scale) * sd(subm_on_toxic$frame$on_toxic_logit_scale)[1] 0.94033613.4 SEM plots
knitr::include_graphics(paste0(getwd(), "/graphs/manuscript1_figures/SEM_simplified.png"), rel_path = FALSE)
4 ISSF parameters
4.1 Selection (immigration) strength state 1
4.1.1 Model summary
m_select_s1_full <- glmmTMB(
s1.less_toxic_end.est ~
cat_pre_wt_log_scale * (beta + var_trt) +
(1|session_id),
data = d3,
); summary(m_select_s1_full) Family: gaussian ( identity )
Formula:
s1.less_toxic_end.est ~ cat_pre_wt_log_scale * (beta + var_trt) +
(1 | session_id)
Data: d3
AIC BIC logLik deviance df.resid
113.7 137.7 -46.9 93.7 71
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 3.344e-11 5.782e-06
Residual 1.862e-01 4.315e-01
Number of obs: 81, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.186
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.035699 0.098240 -0.363 0.716
cat_pre_wt_log_scale 0.124186 0.114435 1.085 0.278
beta0 0.095757 0.119051 0.804 0.421
beta5 -0.036750 0.127318 -0.289 0.773
var_trtlow_var 0.057820 0.100926 0.573 0.567
cat_pre_wt_log_scale:beta0 0.061099 0.137359 0.445 0.656
cat_pre_wt_log_scale:beta5 0.007268 0.146519 0.050 0.960
cat_pre_wt_log_scale:var_trtlow_var -0.078628 0.119191 -0.660 0.509car::Anova(m_select_s1_full, type = "III")Analysis of Deviance Table (Type III Wald chisquare tests)
Response: s1.less_toxic_end.est
Chisq Df Pr(>Chisq)
(Intercept) 0.1320 1 0.7163
cat_pre_wt_log_scale 1.1777 1 0.2778
beta 1.2190 2 0.5436
var_trt 0.3282 1 0.5667
cat_pre_wt_log_scale:beta 0.2245 2 0.8938
cat_pre_wt_log_scale:var_trt 0.4352 1 0.5095m_select_s1 <- glmmTMB(
s1.less_toxic_end.est ~
cat_pre_wt_log_scale + beta + var_trt +
(1|session_id),
data = d3,
); summary(m_select_s1) Family: gaussian ( identity )
Formula:
s1.less_toxic_end.est ~ cat_pre_wt_log_scale + beta + var_trt +
(1 | session_id)
Data: d3
AIC BIC logLik deviance df.resid
108.3 125.1 -47.1 94.3 74
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 3.315e-11 5.758e-06
Residual 1.875e-01 4.331e-01
Number of obs: 81, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.188
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.03455 0.09761 -0.354 0.7234
cat_pre_wt_log_scale 0.10584 0.05788 1.829 0.0675 .
beta0 0.11265 0.11575 0.973 0.3304
beta5 -0.02728 0.11988 -0.228 0.8200
var_trtlow_var 0.03783 0.09662 0.392 0.6954
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(m_select_s1, type = "II")Analysis of Deviance Table (Type II Wald chisquare tests)
Response: s1.less_toxic_end.est
Chisq Df Pr(>Chisq)
cat_pre_wt_log_scale 3.3438 1 0.06746 .
beta 1.5799 2 0.45387
var_trt 0.1533 1 0.69540
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 14.1.2 Post-hoc contrasts
contrast_by_pre_wt(m_select_s1, "cat_pre_wt_log_scale")$emmeans
cat_pre_wt_log_scale emmean SE df lower.CL upper.CL
-2 -0.199 0.139 74 -0.47656 0.0788
2 0.225 0.112 74 0.00197 0.4471
Results are averaged over the levels of: beta, var_trt
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio
(cat_pre_wt_log_scale-2) - cat_pre_wt_log_scale2 -0.423 0.232 74 -1.829
p.value
0.0715
Results are averaged over the levels of: beta, var_trt 4.1.3 Misc stats
r2(m_select_s1, tolerance = 10^-10)Random effect variances not available. Returned R2 does not account for random effects.# R2 for Mixed Models
Conditional R2: NA
Marginal R2: 0.0574.1.4 Check model
check_mod(m_select_s1)
Warning: Maximum value of original data is not included in the
replicated data.
Model may not capture the variation of the data.
4.1.5 Plot
marginal_effects(m_select_s1, terms = c("cat_pre_wt_log_scale")) %>%
ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = exp(yhat))) +
geom_hline(aes(yintercept = 1), color = "tomato", linetype = "dashed", linewidth = 1) +
geom_ribbon(aes(ymax = exp(upper), ymin = exp(lower)), alpha = 0.2) +
geom_line(linewidth = 2) +
geom_point(
data = m_select_s1$frame,
aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = exp(s1.less_toxic_end.est)),
size = 3, alpha = 0.8
) +
theme_bw(base_size = 15) +
theme(legend.position = "top") +
scale_y_continuous(trans = "log10") +
scale_x_continuous(trans = "log10", labels = fancy_scientific) +
labs(x = "Pre-weight (g)", y = "Less toxic diet selection strength (odds ratio)")
4.2 Selection (immigration) strength state 2
4.2.1 Model summary
m_select_s2_full <- glmmTMB(
s2.less_toxic_end.est ~
cat_pre_wt_log_scale * (beta + var_trt) +
(1|session_id),
data = d3 %>%
filter(s2.less_toxic_end.se < 10) # Filter out estimates that did not converge
); summary(m_select_s2_full) Family: gaussian ( identity )
Formula:
s2.less_toxic_end.est ~ cat_pre_wt_log_scale * (beta + var_trt) +
(1 | session_id)
Data: d3 %>% filter(s2.less_toxic_end.se < 10)
AIC BIC logLik deviance df.resid
156.8 180.0 -68.4 136.8 65
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 9.311e-11 9.649e-06
Residual 3.627e-01 6.023e-01
Number of obs: 75, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.363
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.145871 0.138823 1.051 0.293
cat_pre_wt_log_scale 0.076065 0.163210 0.466 0.641
beta0 -0.051651 0.168531 -0.306 0.759
beta5 -0.064804 0.188723 -0.343 0.731
var_trtlow_var -0.006117 0.146957 -0.042 0.967
cat_pre_wt_log_scale:beta0 0.215500 0.194328 1.109 0.267
cat_pre_wt_log_scale:beta5 -0.241814 0.223029 -1.084 0.278
cat_pre_wt_log_scale:var_trtlow_var -0.178505 0.177311 -1.007 0.314car::Anova(m_select_s2_full, type = "III")Analysis of Deviance Table (Type III Wald chisquare tests)
Response: s2.less_toxic_end.est
Chisq Df Pr(>Chisq)
(Intercept) 1.1041 1 0.2934
cat_pre_wt_log_scale 0.2172 1 0.6412
beta 0.1483 2 0.9285
var_trt 0.0017 1 0.9668
cat_pre_wt_log_scale:beta 3.9324 2 0.1400
cat_pre_wt_log_scale:var_trt 1.0135 1 0.3141# Refit
m_select_s2 <- glmmTMB(
s2.less_toxic_end.est ~
cat_pre_wt_log_scale + (var_trt + beta) +
(1|session_id),
data = d3 %>% filter(s2.less_toxic_end.se < 10) # Filter out estimates that did not converge
); summary(m_select_s2) Family: gaussian ( identity )
Formula:
s2.less_toxic_end.est ~ cat_pre_wt_log_scale + (var_trt + beta) +
(1 | session_id)
Data: d3 %>% filter(s2.less_toxic_end.se < 10)
AIC BIC logLik deviance df.resid
154.8 171.0 -70.4 140.8 68
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 6.274e-11 7.921e-06
Residual 3.827e-01 6.186e-01
Number of obs: 75, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.383
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.141025 0.141228 0.999 0.318
cat_pre_wt_log_scale 0.004266 0.086221 0.050 0.961
var_trtlow_var -0.025702 0.144845 -0.177 0.859
beta0 0.005924 0.166898 0.035 0.972
beta5 -0.119943 0.183479 -0.654 0.513car::Anova(m_select_s2, type = "II")Analysis of Deviance Table (Type II Wald chisquare tests)
Response: s2.less_toxic_end.est
Chisq Df Pr(>Chisq)
cat_pre_wt_log_scale 0.0024 1 0.9605
var_trt 0.0315 1 0.8592
beta 0.5632 2 0.75464.2.2 Post-hoc contrasts
emmeans::emmeans(m_select_s2, ~1) 1 emmean SE df lower.CL upper.CL
overall 0.0913 0.0724 68 -0.0532 0.236
Results are averaged over the levels of: var_trt, beta
Confidence level used: 0.95 4.2.3 Misc stats
r2(m_select_s2, tolerance = 10^-10)Random effect variances not available. Returned R2 does not account for random effects.# R2 for Mixed Models
Conditional R2: NA
Marginal R2: 0.0084.2.4 Check model
# Diagnostic plot looks better with m_select_s2 than m_select_s2_2
check_mod(m_select_s2)
Warning: Minimum value of original data is not included in the
replicated data.
Model may not capture the variation of the data.
4.2.5 Plot
marginal_effects(m_select_s2, terms = c("cat_pre_wt_log_scale")) %>%
ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = exp(yhat))) +
geom_hline(aes(yintercept = 1), color = "tomato", linetype = "dashed", linewidth = 1) +
geom_ribbon(aes(ymax = exp(upper), ymin = exp(lower)), alpha = 0.2) +
geom_line(linewidth = 2) +
geom_point(
data = m_select_s2$frame,
aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = exp(s2.less_toxic_end.est)),
size = 3, alpha = 0.8
) +
theme_bw(base_size = 15) +
theme(legend.position = "top") +
scale_y_continuous(trans = "log10") +
scale_x_continuous(trans = "log10", labels = fancy_scientific) +
labs(x = "Pre-weight (g)", y = "Less toxic diet selection strength (odds ratio)")
4.3 Arrestment state 1
4.3.1 Model summary
m_arrest_s1_full <- glmmTMB(
log(k1) ~
cat_pre_wt_log_scale * (beta + var_trt) +
(1|session_id),
data = d3,
); summary(m_arrest_s1_full) Family: gaussian ( identity )
Formula:
log(k1) ~ cat_pre_wt_log_scale * (beta + var_trt) + (1 | session_id)
Data: d3
AIC BIC logLik deviance df.resid
202.2 226.2 -91.1 182.2 71
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 0.000537 0.02317
Residual 0.554954 0.74495
Number of obs: 81, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.555
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.08923 0.17091 -0.522 0.6016
cat_pre_wt_log_scale 0.04337 0.20406 0.212 0.8317
beta0 0.08266 0.20567 0.402 0.6877
beta5 -0.09266 0.21995 -0.421 0.6736
var_trtlow_var -0.04815 0.17442 -0.276 0.7825
cat_pre_wt_log_scale:beta0 0.10559 0.23857 0.443 0.6581
cat_pre_wt_log_scale:beta5 0.63298 0.25306 2.501 0.0124 *
cat_pre_wt_log_scale:var_trtlow_var 0.09904 0.20584 0.481 0.6304
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(m_arrest_s1_full, type = "III")Analysis of Deviance Table (Type III Wald chisquare tests)
Response: log(k1)
Chisq Df Pr(>Chisq)
(Intercept) 0.2726 1 0.60161
cat_pre_wt_log_scale 0.0452 1 0.83170
beta 0.6378 2 0.72694
var_trt 0.0762 1 0.78251
cat_pre_wt_log_scale:beta 6.7997 2 0.03338 *
cat_pre_wt_log_scale:var_trt 0.2315 1 0.63039
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m_arrest_s1 <- glmmTMB(
log(k1) ~
cat_pre_wt_log_scale * beta + var_trt +
(1|session_id),
data = d3,
); summary(m_arrest_s1) Family: gaussian ( identity )
Formula:
log(k1) ~ cat_pre_wt_log_scale * beta + var_trt + (1 | session_id)
Data: d3
AIC BIC logLik deviance df.resid
200.5 222.0 -91.2 182.5 72
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 0.0001243 0.01115
Residual 0.5569498 0.74629
Number of obs: 81, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.557
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.09552 0.17047 -0.560 0.5753
cat_pre_wt_log_scale 0.09510 0.17338 0.548 0.5834
beta0 0.07431 0.20529 0.362 0.7174
beta5 -0.10004 0.21979 -0.455 0.6490
var_trtlow_var -0.02513 0.16800 -0.150 0.8811
cat_pre_wt_log_scale:beta0 0.11796 0.23764 0.496 0.6196
cat_pre_wt_log_scale:beta5 0.61691 0.25131 2.455 0.0141 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(m_arrest_s1, type = "II")Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(k1)
Chisq Df Pr(>Chisq)
cat_pre_wt_log_scale 6.8531 1 0.008849 **
beta 0.3005 2 0.860488
var_trt 0.0224 1 0.881103
cat_pre_wt_log_scale:beta 6.5451 2 0.037911 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(m_arrest_s1, type = "III")Analysis of Deviance Table (Type III Wald chisquare tests)
Response: log(k1)
Chisq Df Pr(>Chisq)
(Intercept) 0.3139 1 0.57527
cat_pre_wt_log_scale 0.3008 1 0.58336
beta 0.6279 2 0.73057
var_trt 0.0224 1 0.88110
cat_pre_wt_log_scale:beta 6.5451 2 0.03791 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 14.3.2 Post-hoc contrasts
contrast_by_pre_wt(m_arrest_s1, "beta")$emmeans
cat_pre_wt_log_scale = -2:
beta emmean SE df lower.CL upper.CL
-5 -0.2983 0.409 72 -1.114 0.517
0 -0.4599 0.432 72 -1.320 0.401
5 -1.6321 0.492 72 -2.612 -0.652
cat_pre_wt_log_scale = 2:
beta emmean SE df lower.CL upper.CL
-5 0.0821 0.340 72 -0.597 0.761
0 0.3923 0.370 72 -0.345 1.129
5 1.2159 0.356 72 0.505 1.926
Results are averaged over the levels of: var_trt
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
cat_pre_wt_log_scale = -2:
contrast estimate SE df t.ratio p.value
(beta-5) - beta0 0.162 0.559 72 0.289 0.9550
(beta-5) - beta5 1.334 0.613 72 2.178 0.0819
beta0 - beta5 1.172 0.615 72 1.905 0.1446
cat_pre_wt_log_scale = 2:
contrast estimate SE df t.ratio p.value
(beta-5) - beta0 -0.310 0.473 72 -0.656 0.7894
(beta-5) - beta5 -1.134 0.476 72 -2.381 0.0514
beta0 - beta5 -0.824 0.483 72 -1.705 0.2104
Results are averaged over the levels of: var_trt
Results are given on the log (not the response) scale.
P value adjustment: tukey method for comparing a family of 3 estimates 4.3.3 Misc stats
r2(m_arrest_s1, tolerance = 10^-10)# R2 for Mixed Models
Conditional R2: 0.175
Marginal R2: 0.1754.3.4 Check model
check_mod(m_arrest_s1)
Warning: Maximum value of original data is not included in the
replicated data.
Model may not capture the variation of the data.
4.3.5 Plot
marginal_effects(m_arrest_s1, terms = c("cat_pre_wt_log_scale", "beta")) %>%
ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = (yhat), color = beta)) +
geom_ribbon(aes(ymax = (upper), ymin = lower, fill = beta, color = NULL), alpha = 0.2) +
geom_line(linewidth = 2) +
geom_point(
data = m_arrest_s1$frame,
aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = exp(`log(k1)`)),
size = 3, alpha = 0.5
) +
theme_bw(base_size = 15) +
geom_hline(aes(yintercept = 1), color = "black", linetype = "dashed", linewidth = 1) +
theme(legend.position = "top") +
scale_y_continuous(trans = "log10") +
scale_x_continuous(trans = "log10", labels = fancy_scientific) +
scale_color_brewer(type = "qual", aesthetics = c("color","fill")) +
labs(x = "Pre-weight (g)", y = "Less toxic diet arresetment strength (ratio)",
color = expression(beta), fill = expression(beta))
4.4 Arrestment state 2
4.4.1 Model summary
m_arrest_s2_full <- glmmTMB(
log(k2) ~
cat_pre_wt_log_scale * (beta + var_trt) +
(1|session_id),
data = d3,
); summary(m_arrest_s2_full) Family: gaussian ( identity )
Formula:
log(k2) ~ cat_pre_wt_log_scale * (beta + var_trt) + (1 | session_id)
Data: d3
AIC BIC logLik deviance df.resid
154 178 -67 134 71
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 0.05334 0.2310
Residual 0.28154 0.5306
Number of obs: 81, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.282
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.20619 0.15939 1.294 0.1958
cat_pre_wt_log_scale -0.13519 0.18141 -0.745 0.4561
beta0 -0.09484 0.14743 -0.643 0.5201
beta5 -0.18572 0.15717 -1.182 0.2373
var_trtlow_var -0.25005 0.12598 -1.985 0.0472 *
cat_pre_wt_log_scale:beta0 0.07897 0.16980 0.465 0.6419
cat_pre_wt_log_scale:beta5 0.10713 0.18053 0.593 0.5529
cat_pre_wt_log_scale:var_trtlow_var 0.09679 0.14916 0.649 0.5164
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(m_arrest_s2_full, type = "III")Analysis of Deviance Table (Type III Wald chisquare tests)
Response: log(k2)
Chisq Df Pr(>Chisq)
(Intercept) 1.6735 1 0.19579
cat_pre_wt_log_scale 0.5554 1 0.45612
beta 1.4088 2 0.49442
var_trt 3.9397 1 0.04716 *
cat_pre_wt_log_scale:beta 0.4036 2 0.81725
cat_pre_wt_log_scale:var_trt 0.4211 1 0.51637
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m_arrest_s2 <- glmmTMB(
log(k2) ~
cat_pre_wt_log_scale + beta + var_trt +
(1|session_id),
data = d3,
); summary(m_arrest_s2) Family: gaussian ( identity )
Formula:
log(k2) ~ cat_pre_wt_log_scale + beta + var_trt + (1 | session_id)
Data: d3
AIC BIC logLik deviance df.resid
148.8 165.6 -67.4 134.8 74
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 0.04395 0.2096
Residual 0.28700 0.5357
Number of obs: 81, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.287
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.19258 0.15345 1.255 0.2095
cat_pre_wt_log_scale -0.01354 0.13190 -0.103 0.9182
beta0 -0.08219 0.14427 -0.570 0.5689
beta5 -0.16897 0.14871 -1.136 0.2558
var_trtlow_var -0.23314 0.12060 -1.933 0.0532 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(m_arrest_s2, type = "II")Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(k2)
Chisq Df Pr(>Chisq)
cat_pre_wt_log_scale 0.0105 1 0.91823
beta 1.2923 2 0.52406
var_trt 3.7367 1 0.05323 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 14.4.2 Post-hoc contrasts
emmeans::emmeans(m_arrest_s2, ~ "var_trt") var_trt emmean SE df lower.CL upper.CL
high_var 0.105 0.132 74 -0.158 0.369
low_var -0.128 0.137 74 -0.400 0.145
Results are averaged over the levels of: beta
Results are given on the log (not the response) scale.
Confidence level used: 0.95 4.4.3 Misc stats
r2(m_arrest_s2, tolerance = 10^-10)# R2 for Mixed Models
Conditional R2: 0.176
Marginal R2: 0.0504.4.4 Check model
check_mod(m_arrest_s2)
Warning: Minimum value of original data is not included in the
replicated data.
Model may not capture the variation of the data.
4.4.5 Plot
marginal_effects(m_arrest_s2, terms = c("cat_pre_wt_log_scale", "beta")) %>%
ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = (yhat), color = beta)) +
geom_ribbon(aes(ymax = (upper), ymin = lower, fill = beta, color = NULL), alpha = 0.2) +
geom_line(linewidth = 2) +
geom_point(
data = m_arrest_s2$frame,
aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = exp(`log(k2)`)),
size = 3, alpha = 0.5
) +
theme_bw(base_size = 15) +
geom_hline(aes(yintercept = 1), color = "black", linetype = "dashed", linewidth = 1) +
theme(legend.position = "top") +
scale_y_continuous(trans = "log10") +
scale_x_continuous(trans = "log10", labels = fancy_scientific) +
scale_color_brewer(type = "qual", aesthetics = c("color","fill")) +
labs(x = "Pre-weight (g)", y = "Less toxic diet arresetment strength (ratio)",
color = expression(beta), fill = expression(beta))
4.5 Toxic diet step length state 1
4.5.1 Model summary
m_sl_pred_s1_full <- glmmTMB(
log(sl_mean_pred1) ~
cat_pre_wt_log_scale * (beta + var_trt) +
(1|session_id),
data = d3,
); summary(m_sl_pred_s1_full) Family: gaussian ( identity )
Formula:
log(sl_mean_pred1) ~ cat_pre_wt_log_scale * (beta + var_trt) +
(1 | session_id)
Data: d3
AIC BIC logLik deviance df.resid
262.6 286.5 -121.3 242.6 71
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 8.455e-10 2.908e-05
Residual 1.170e+00 1.082e+00
Number of obs: 81, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 1.17
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.98102 0.24627 16.165 <2e-16 ***
cat_pre_wt_log_scale -0.03586 0.28686 -0.125 0.901
beta0 0.19843 0.29844 0.665 0.506
beta5 0.43478 0.31916 1.362 0.173
var_trtlow_var 0.34703 0.25300 1.372 0.170
cat_pre_wt_log_scale:beta0 0.10587 0.34433 0.307 0.758
cat_pre_wt_log_scale:beta5 -0.12070 0.36729 -0.329 0.742
cat_pre_wt_log_scale:var_trtlow_var 0.22539 0.29879 0.754 0.451
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(m_sl_pred_s1_full, type = "III")Analysis of Deviance Table (Type III Wald chisquare tests)
Response: log(sl_mean_pred1)
Chisq Df Pr(>Chisq)
(Intercept) 261.3177 1 <2e-16 ***
cat_pre_wt_log_scale 0.0156 1 0.9005
beta 1.8565 2 0.3952
var_trt 1.8814 1 0.1702
cat_pre_wt_log_scale:beta 0.3619 2 0.8345
cat_pre_wt_log_scale:var_trt 0.5690 1 0.4506
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m_sl_pred_s1 <- glmmTMB(
log(sl_mean_pred1) ~
cat_pre_wt_log_scale + beta + var_trt +
(1|session_id),
data = d3,
); summary(m_sl_pred_s1) Family: gaussian ( identity )
Formula: log(sl_mean_pred1) ~ cat_pre_wt_log_scale + beta + var_trt +
(1 | session_id)
Data: d3
AIC BIC logLik deviance df.resid
257.8 274.6 -121.9 243.8 74
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 7.726e-10 2.78e-05
Residual 1.188e+00 1.09e+00
Number of obs: 81, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 1.19
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.95860 0.24562 16.117 <2e-16 ***
cat_pre_wt_log_scale 0.08645 0.14566 0.594 0.5528
beta0 0.20629 0.29126 0.708 0.4788
beta5 0.36034 0.30166 1.195 0.2323
var_trtlow_var 0.41258 0.24315 1.697 0.0897 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(m_sl_pred_s1, type = "II")Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(sl_mean_pred1)
Chisq Df Pr(>Chisq)
cat_pre_wt_log_scale 0.3523 1 0.55282
beta 1.4497 2 0.48440
var_trt 2.8793 1 0.08973 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 14.5.2 Post-hoc contrasts
emmeans::emmeans(m_sl_pred_s1, ~ "var_trt") var_trt emmean SE df lower.CL upper.CL
high_var 4.17 0.172 74 3.83 4.51
low_var 4.58 0.171 74 4.24 4.92
Results are averaged over the levels of: beta
Results are given on the log (not the response) scale.
Confidence level used: 0.95 emmeans::emmeans(m_sl_pred_s1, ~ 1) 1 emmean SE df lower.CL upper.CL
overall 4.38 0.121 74 4.13 4.62
Results are averaged over the levels of: beta, var_trt
Results are given on the log (not the response) scale.
Confidence level used: 0.95 4.5.3 Misc stats
r2(m_sl_pred_s1, tolerance = 10^-10)# R2 for Mixed Models
Conditional R2: 0.053
Marginal R2: 0.0534.5.4 Check model
check_mod(m_sl_pred_s1)
4.5.5 Plot
marginal_effects(m_sl_pred_s1, terms = c("cat_pre_wt_log_scale")) %>%
ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = (yhat) / 1000 * 12)) +
geom_ribbon(aes(ymax = (upper)/ 1000 * 12, ymin = lower/ 1000 * 12), alpha = 0.2) +
geom_line(linewidth = 2) +
geom_point(
data = m_sl_pred_s1$frame,
aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale),
y = exp(`log(sl_mean_pred1)`) / 1000 * 12),
size = 3, alpha = 0.5
) +
theme_bw(base_size = 15) +
theme(legend.position = "top") +
scale_y_continuous(trans = "log10", labels = fancy_scientific) +
scale_x_continuous(trans = "log10", labels = fancy_scientific) +
scale_color_brewer(type = "qual", aesthetics = c("color","fill")) +
labs(x = "Cat pre-weight (g)", y = "Selection free step length \non toxic diet (cm)")
4.6 Toxic diet step length state 2
4.6.1 Model summary
m_sl_pred_s2_full <- glmmTMB(
log(sl_mean_pred2) ~
cat_pre_wt_log_scale * (beta + var_trt) +
(1|session_id),
data = d3,
); summary(m_sl_pred_s2_full) Family: gaussian ( identity )
Formula:
log(sl_mean_pred2) ~ cat_pre_wt_log_scale * (beta + var_trt) +
(1 | session_id)
Data: d3
AIC BIC logLik deviance df.resid
71.5 95.5 -25.8 51.5 71
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 6.870e-12 2.621e-06
Residual 1.106e-01 3.326e-01
Number of obs: 81, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.111
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.014673 0.075709 26.611 <2e-16 ***
cat_pre_wt_log_scale 0.036250 0.088189 0.411 0.6810
beta0 -0.159917 0.091747 -1.743 0.0813 .
beta5 -0.083506 0.098117 -0.851 0.3947
var_trtlow_var -0.068098 0.077779 -0.876 0.3813
cat_pre_wt_log_scale:beta0 0.156273 0.105856 1.476 0.1399
cat_pre_wt_log_scale:beta5 -0.131417 0.112915 -1.164 0.2445
cat_pre_wt_log_scale:var_trtlow_var -0.008615 0.091854 -0.094 0.9253
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Somehow the Wald chisquare test is significant.
car::Anova(m_sl_pred_s2_full, type = "III")Analysis of Deviance Table (Type III Wald chisquare tests)
Response: log(sl_mean_pred2)
Chisq Df Pr(>Chisq)
(Intercept) 708.1363 1 < 2e-16 ***
cat_pre_wt_log_scale 0.1690 1 0.68103
beta 3.0394 2 0.21878
var_trt 0.7666 1 0.38128
cat_pre_wt_log_scale:beta 6.2659 2 0.04359 *
cat_pre_wt_log_scale:var_trt 0.0088 1 0.92528
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m_sl_pred_s2 <- glmmTMB(
log(sl_mean_pred2) ~
cat_pre_wt_log_scale * beta + var_trt +
(1|session_id),
data = d3,
); summary(m_sl_pred_s2) Family: gaussian ( identity )
Formula: log(sl_mean_pred2) ~ cat_pre_wt_log_scale * beta + var_trt +
(1 | session_id)
Data: d3
AIC BIC logLik deviance df.resid
69.5 91.1 -25.8 51.5 72
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
session_id (Intercept) 7.862e-12 2.804e-06
Residual 1.106e-01 3.326e-01
Number of obs: 81, groups: session_id, 5
Dispersion estimate for gaussian family (sigma^2): 0.111
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.01525 0.07547 26.704 <2e-16 ***
cat_pre_wt_log_scale 0.03168 0.07348 0.431 0.6664
beta0 -0.15920 0.09143 -1.741 0.0817 .
beta5 -0.08285 0.09787 -0.847 0.3973
var_trtlow_var -0.07011 0.07476 -0.938 0.3483
cat_pre_wt_log_scale:beta0 0.15516 0.10519 1.475 0.1402
cat_pre_wt_log_scale:beta5 -0.13003 0.11195 -1.162 0.2454
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(m_sl_pred_s2, type = "III")Analysis of Deviance Table (Type III Wald chisquare tests)
Response: log(sl_mean_pred2)
Chisq Df Pr(>Chisq)
(Intercept) 713.1060 1 < 2e-16 ***
cat_pre_wt_log_scale 0.1858 1 0.66641
beta 3.0330 2 0.21948
var_trt 0.8796 1 0.34830
cat_pre_wt_log_scale:beta 6.4980 2 0.03881 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 14.6.2 Post-hoc contrasts
contrast_by_pre_wt(m_sl_pred_s2, "beta")$emmeans
cat_pre_wt_log_scale = -2:
beta emmean SE df lower.CL upper.CL
-5 1.92 0.174 72 1.57 2.26
0 1.45 0.177 72 1.09 1.80
5 2.09 0.209 72 1.68 2.51
cat_pre_wt_log_scale = 2:
beta emmean SE df lower.CL upper.CL
-5 2.04 0.146 72 1.75 2.33
0 2.19 0.150 72 1.90 2.49
5 1.70 0.153 72 1.40 2.00
Results are averaged over the levels of: var_trt
Results are given on the log (not the response) scale.
Confidence level used: 0.95
$contrasts
cat_pre_wt_log_scale = -2:
contrast estimate SE df t.ratio p.value
(beta-5) - beta0 0.470 0.248 72 1.891 0.1486
(beta-5) - beta5 -0.177 0.273 72 -0.650 0.7932
beta0 - beta5 -0.647 0.274 72 -2.360 0.0540
cat_pre_wt_log_scale = 2:
contrast estimate SE df t.ratio p.value
(beta-5) - beta0 -0.151 0.209 72 -0.724 0.7503
(beta-5) - beta5 0.343 0.212 72 1.616 0.2454
beta0 - beta5 0.494 0.214 72 2.310 0.0608
Results are averaged over the levels of: var_trt
Results are given on the log (not the response) scale.
P value adjustment: tukey method for comparing a family of 3 estimates emmeans::emmeans(m_sl_pred_s2, ~ 1) 1 emmean SE df lower.CL upper.CL
overall 1.91 0.0372 72 1.84 1.98
Results are averaged over the levels of: beta, var_trt
Results are given on the log (not the response) scale.
Confidence level used: 0.95 4.6.3 Misc stats
r2(m_sl_pred_s2, tolerance = 10^-10)Random effect variances not available. Returned R2 does not account for random effects.# R2 for Mixed Models
Conditional R2: NA
Marginal R2: 0.1204.6.4 Check model
check_mod(m_sl_pred_s2)
4.6.5 Plot
marginal_effects(m_sl_pred_s2, terms = c("cat_pre_wt_log_scale", "beta")) %>%
ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale),
y = (yhat) / 1000 * 12,
color = beta)) +
geom_ribbon(aes(ymax = (upper)/ 1000 * 12, ymin = lower/ 1000 * 12,
fill = beta, color = NULL),
alpha = 0.2) +
geom_line(linewidth = 2) +
geom_point(
data = m_sl_pred_s2$frame,
aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale),
y = exp(`log(sl_mean_pred2)`) / 1000 * 12),
size = 3, alpha = 0.5
) +
theme_bw(base_size = 15) +
theme(legend.position = "top") +
scale_y_continuous(trans = "log10") +
scale_x_continuous(trans = "log10", labels = fancy_scientific) +
scale_color_brewer(type = "qual", aesthetics = c("color","fill")) +
labs(x = "Pre-weight (g)", y = "Selection free step length \non toxic diet (cm)")
5 ISSF simulation
5.1 Observed means
5.1.1 Summary
# m <- brm(
# on_toxic ~
# beta * cat_pre_wt_log_scale + var_trt +
# (1|session_id),
# family = Beta(),
# data = d %>%
# filter(var_trt != "constant"),
# chains = 4,
# cores = 4,
# prior = c(
# set_prior("normal(0,1.5)", class = "b"),
# set_prior("normal(0,1)", class = "Intercept"),
# set_prior("cauchy(0,0.5)", class = "sd"),
# set_prior("gamma(0.01, 0.01)", class = "phi")
# ),
# iter = 4000,
# control = list(adapt_delta = 0.95)
# )
# saveRDS(m, "invisible/fitted_models/on_toxic_brm.rds")
m <- readRDS("invisible/fitted_models/on_toxic_brm.rds")
m Family: beta
Links: mu = logit; phi = identity
Formula: on_toxic ~ beta * cat_pre_wt_log_scale + var_trt + (1 | session_id)
Data: d %>% filter(var_trt != "constant") (Number of observations: 97)
Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
total post-warmup draws = 8000
Group-Level Effects:
~session_id (Number of levels: 5)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.24 0.19 0.01 0.72 1.00 2231 3150
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept -0.73 0.23 -1.19 -0.27 1.00 4467
beta0 0.14 0.22 -0.29 0.56 1.00 7373
beta5 -0.20 0.24 -0.66 0.28 1.00 6607
cat_pre_wt_log_scale 0.10 0.20 -0.28 0.49 1.00 4581
var_trtlow_var 0.36 0.18 0.01 0.72 1.00 8461
beta0:cat_pre_wt_log_scale -0.30 0.24 -0.76 0.18 1.00 5456
beta5:cat_pre_wt_log_scale -0.72 0.29 -1.29 -0.16 1.00 5511
Tail_ESS
Intercept 4001
beta0 5574
beta5 5527
cat_pre_wt_log_scale 5668
var_trtlow_var 5568
beta0:cat_pre_wt_log_scale 5516
beta5:cat_pre_wt_log_scale 6012
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
phi 4.72 0.66 3.51 6.09 1.00 6699 5717
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).5.1.2 Check model
brms:::plot.brmsfit(m)
brms::pp_check(m,ndraws = 100)
5.1.3 Plot
brms::conditional_effects(m,effects = c("cat_pre_wt_log_scale:beta"), re_formula = NA)[[1]] %>%
ggplot(aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = estimate__, color = beta)) +
geom_ribbon(aes(ymax = upper__, ymin = lower__, color = NULL, fill = beta), alpha = 0.2) +
geom_line(size = 2) +
geom_point(
data = m$data,
aes(x = unscalelog(d$cat_pre_wt_log)(cat_pre_wt_log_scale), y = on_toxic),
size = 5,
alpha = 0.5
) +
theme_bw(base_size = 15) +
scale_x_continuous(trans = "log10", labels = fancy_scientific) +
scale_color_brewer(type = "qual", aesthetics = c("fill", "color")) +
labs(x = "Pre-weight (g)", y = "Proprtion of time on toxic diet", color = expression(beta), fill = expression(beta))
5.1.4 Observed effect sizes
# Compute observed effect size
epred_draws(m, expand.grid("cat_pre_wt_log_scale" = c(-2, 2),
"beta" = c(-5, 5),
"var_trt" = c("low_var","high_var"),
"session_id" = "foo")) %>%
group_by(cat_pre_wt_log_scale, beta, draw) %>%
summarise(val = mean(val)) %>%
group_by(
cat_pre_wt_log_scale, draw
) %>%
summarise(
delta = diff(val)
) %>%
group_by(cat_pre_wt_log_scale) %>%
do(as.data.frame(t(summarise_vec(.$delta))))# A tibble: 2 × 6
# Groups: cat_pre_wt_log_scale [2]
cat_pre_wt_log_scale mean median sd lower upper
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 -2 0.280 0.284 0.151 -0.0251 0.563
2 2 -0.285 -0.284 0.0970 -0.479 -0.09915.2 Group means
out <- read_csv("simulation/move_rules_sim.csv")
out <- out %>%
mutate(
scale = round(scale / 1000 * 12, 2),
rss = as.factor(round(exp(rss), 2))
)
out %>%
group_by(b, rss, scale, k) %>%
summarise(
mean = mean(on_toxic), se = se(on_toxic), sd = sd(on_toxic),
lower = quantile(on_toxic, probs = 0.025), upper = quantile(on_toxic, probs = 0.975)
) %>%
mutate(
z = (mean - 0.5) / se,
P = pnorm(abs(z), lower.tail = FALSE) * 2
) %>%
knitr::kable()| b | rss | scale | k | mean | se | sd | lower | upper | z | P |
|---|---|---|---|---|---|---|---|---|---|---|
| -5 | 0.8 | 0.1 | 0.25 | 0.9039186 | 0.0018517 | 0.0370341 | 0.8151349 | 0.9671079 | 218.1331185 | 0.0000000 |
| -5 | 0.8 | 0.1 | 1.00 | 0.5905470 | 0.0069139 | 0.1382788 | 0.3143606 | 0.8511489 | 13.0962854 | 0.0000000 |
| -5 | 0.8 | 0.1 | 4.00 | 0.1331119 | 0.0052640 | 0.1052799 | 0.0029970 | 0.3729520 | -69.6976278 | 0.0000000 |
| -5 | 0.8 | 0.3 | 0.25 | 0.8059266 | 0.0017740 | 0.0354790 | 0.7432567 | 0.8761239 | 172.4549517 | 0.0000000 |
| -5 | 0.8 | 0.3 | 1.00 | 0.5962987 | 0.0038047 | 0.0760946 | 0.4274476 | 0.7482517 | 25.3102497 | 0.0000000 |
| -5 | 0.8 | 0.3 | 4.00 | 0.1527298 | 0.0033212 | 0.0664248 | 0.0459041 | 0.3127872 | -104.5604328 | 0.0000000 |
| -5 | 0.8 | 1.0 | 0.25 | 0.6634316 | 0.0011939 | 0.0238778 | 0.6153596 | 0.7083666 | 136.8900920 | 0.0000000 |
| -5 | 0.8 | 1.0 | 1.00 | 0.5805495 | 0.0016245 | 0.0324908 | 0.5124126 | 0.6473776 | 49.5828564 | 0.0000000 |
| -5 | 0.8 | 1.0 | 4.00 | 0.2898601 | 0.0024510 | 0.0490205 | 0.1928072 | 0.3846653 | -85.7354275 | 0.0000000 |
| -5 | 1 | 0.1 | 0.25 | 0.8788961 | 0.0023185 | 0.0463708 | 0.7831918 | 0.9590410 | 163.4202417 | 0.0000000 |
| -5 | 1 | 0.1 | 1.00 | 0.4908342 | 0.0061322 | 0.1226433 | 0.2497253 | 0.7343656 | -1.4947137 | 0.1349892 |
| -5 | 1 | 0.1 | 4.00 | 0.0950150 | 0.0043984 | 0.0879685 | 0.0009990 | 0.3289960 | -92.0750307 | 0.0000000 |
| -5 | 1 | 0.3 | 0.25 | 0.7522952 | 0.0020909 | 0.0418177 | 0.6762737 | 0.8292208 | 120.6642601 | 0.0000000 |
| -5 | 1 | 0.3 | 1.00 | 0.5039910 | 0.0037582 | 0.0751641 | 0.3644605 | 0.6543457 | 1.0619450 | 0.2882606 |
| -5 | 1 | 0.3 | 4.00 | 0.1090909 | 0.0024322 | 0.0486447 | 0.0359141 | 0.2138611 | -160.7200891 | 0.0000000 |
| -5 | 1 | 1.0 | 0.25 | 0.5968731 | 0.0012855 | 0.0257093 | 0.5444306 | 0.6463786 | 75.3604269 | 0.0000000 |
| -5 | 1 | 1.0 | 1.00 | 0.4997977 | 0.0017462 | 0.0349246 | 0.4315185 | 0.5704296 | -0.1158482 | 0.9077728 |
| -5 | 1 | 1.0 | 4.00 | 0.2229146 | 0.0019965 | 0.0399302 | 0.1458541 | 0.2987013 | -138.7848297 | 0.0000000 |
| -5 | 1.25 | 0.1 | 0.25 | 0.8291958 | 0.0029629 | 0.0592574 | 0.7012987 | 0.9380869 | 111.1070747 | 0.0000000 |
| -5 | 1.25 | 0.1 | 1.00 | 0.4059590 | 0.0067969 | 0.1359388 | 0.1637862 | 0.6994755 | -13.8357759 | 0.0000000 |
| -5 | 1.25 | 0.1 | 4.00 | 0.0732118 | 0.0038774 | 0.0775487 | 0.0009990 | 0.2867383 | -110.0696915 | 0.0000000 |
| -5 | 1.25 | 0.3 | 0.25 | 0.6839660 | 0.0022279 | 0.0445578 | 0.6013237 | 0.7672827 | 82.5740974 | 0.0000000 |
| -5 | 1.25 | 0.3 | 1.00 | 0.3934191 | 0.0035980 | 0.0719601 | 0.2627123 | 0.5456044 | -29.6222231 | 0.0000000 |
| -5 | 1.25 | 0.3 | 4.00 | 0.0847577 | 0.0019719 | 0.0394377 | 0.0229520 | 0.1689061 | -210.5814170 | 0.0000000 |
| -5 | 1.25 | 1.0 | 0.25 | 0.5268457 | 0.0013041 | 0.0260821 | 0.4764985 | 0.5764735 | 20.5854906 | 0.0000000 |
| -5 | 1.25 | 1.0 | 1.00 | 0.4200100 | 0.0017184 | 0.0343685 | 0.3496503 | 0.4895854 | -46.5485105 | 0.0000000 |
| -5 | 1.25 | 1.0 | 4.00 | 0.1670554 | 0.0017809 | 0.0356186 | 0.0978272 | 0.2427822 | -186.9499930 | 0.0000000 |
| 0 | 0.8 | 0.1 | 0.25 | 0.9154396 | 0.0020826 | 0.0416514 | 0.8289960 | 0.9800450 | 199.4841251 | 0.0000000 |
| 0 | 0.8 | 0.1 | 1.00 | 0.6122827 | 0.0083901 | 0.1678026 | 0.2497003 | 0.8991009 | 13.3827165 | 0.0000000 |
| 0 | 0.8 | 0.1 | 4.00 | 0.1549850 | 0.0080624 | 0.1612470 | 0.0009990 | 0.5434815 | -42.7933403 | 0.0000000 |
| 0 | 0.8 | 0.3 | 0.25 | 0.8413237 | 0.0018827 | 0.0376546 | 0.7662338 | 0.9171329 | 181.2919850 | 0.0000000 |
| 0 | 0.8 | 0.3 | 1.00 | 0.6001873 | 0.0047834 | 0.0956689 | 0.4034965 | 0.7932817 | 20.9446032 | 0.0000000 |
| 0 | 0.8 | 0.3 | 4.00 | 0.1455420 | 0.0043286 | 0.0865713 | 0.0268981 | 0.3247502 | -81.8881337 | 0.0000000 |
| 0 | 0.8 | 1.0 | 0.25 | 0.6926149 | 0.0014591 | 0.0291815 | 0.6363636 | 0.7472527 | 132.0114712 | 0.0000000 |
| 0 | 0.8 | 1.0 | 1.00 | 0.5832243 | 0.0022441 | 0.0448828 | 0.4984016 | 0.6674575 | 37.0851532 | 0.0000000 |
| 0 | 0.8 | 1.0 | 4.00 | 0.2501199 | 0.0030531 | 0.0610614 | 0.1368631 | 0.3816434 | -81.8455512 | 0.0000000 |
| 0 | 1 | 0.1 | 0.25 | 0.8816234 | 0.0030176 | 0.0603529 | 0.7442308 | 0.9780220 | 126.4640036 | 0.0000000 |
| 0 | 1 | 0.1 | 1.00 | 0.5025924 | 0.0085678 | 0.1713558 | 0.1747752 | 0.8043706 | 0.3025761 | 0.7622129 |
| 0 | 1 | 0.1 | 4.00 | 0.1114461 | 0.0065858 | 0.1317168 | 0.0009990 | 0.4890360 | -58.9983792 | 0.0000000 |
| 0 | 1 | 0.3 | 0.25 | 0.7911538 | 0.0024125 | 0.0482507 | 0.6923077 | 0.8841658 | 120.6837732 | 0.0000000 |
| 0 | 1 | 0.3 | 1.00 | 0.5023701 | 0.0050301 | 0.1006024 | 0.3096653 | 0.6893357 | 0.4711874 | 0.6375069 |
| 0 | 1 | 0.3 | 4.00 | 0.1084366 | 0.0030472 | 0.0609447 | 0.0169081 | 0.2468282 | -128.4979711 | 0.0000000 |
| 0 | 1 | 1.0 | 0.25 | 0.6282917 | 0.0016147 | 0.0322930 | 0.5664086 | 0.6933067 | 79.4547584 | 0.0000000 |
| 0 | 1 | 1.0 | 1.00 | 0.5003621 | 0.0021468 | 0.0429356 | 0.4135864 | 0.5884116 | 0.1686888 | 0.8660415 |
| 0 | 1 | 1.0 | 4.00 | 0.1918981 | 0.0024640 | 0.0492807 | 0.0948801 | 0.2918082 | -125.0394879 | 0.0000000 |
| 0 | 1.25 | 0.1 | 0.25 | 0.8398027 | 0.0036314 | 0.0726281 | 0.6822927 | 0.9620629 | 93.5733645 | 0.0000000 |
| 0 | 1.25 | 0.1 | 1.00 | 0.4036014 | 0.0086244 | 0.1724887 | 0.1177073 | 0.7423576 | -11.1773828 | 0.0000000 |
| 0 | 1.25 | 0.1 | 4.00 | 0.0906968 | 0.0053304 | 0.1066081 | 0.0009990 | 0.4089411 | -76.7865320 | 0.0000000 |
| 0 | 1.25 | 0.3 | 0.25 | 0.7233292 | 0.0028702 | 0.0574040 | 0.6082667 | 0.8372128 | 77.8096042 | 0.0000000 |
| 0 | 1.25 | 0.3 | 1.00 | 0.4058666 | 0.0048732 | 0.0974646 | 0.2227023 | 0.5925325 | -19.3164192 | 0.0000000 |
| 0 | 1.25 | 0.3 | 4.00 | 0.0789535 | 0.0023950 | 0.0479005 | 0.0129620 | 0.2027972 | -175.8004437 | 0.0000000 |
| 0 | 1.25 | 1.0 | 0.25 | 0.5525574 | 0.0017095 | 0.0341894 | 0.4844156 | 0.6173826 | 30.7449156 | 0.0000000 |
| 0 | 1.25 | 1.0 | 1.00 | 0.4180020 | 0.0021056 | 0.0421124 | 0.3356643 | 0.4956543 | -38.9424151 | 0.0000000 |
| 0 | 1.25 | 1.0 | 4.00 | 0.1420055 | 0.0020180 | 0.0403607 | 0.0629121 | 0.2238262 | -177.3975346 | 0.0000000 |
| 5 | 0.8 | 0.1 | 0.25 | 0.9130170 | 0.0044098 | 0.0881967 | 0.6540210 | 0.9990010 | 93.6581027 | 0.0000000 |
| 5 | 0.8 | 0.1 | 1.00 | 0.5673002 | 0.0143028 | 0.2860567 | 0.0029471 | 1.0000000 | 4.7053742 | 0.0000025 |
| 5 | 0.8 | 0.1 | 4.00 | 0.2545105 | 0.0151281 | 0.3025619 | 0.0009990 | 1.0000000 | -16.2273900 | 0.0000000 |
| 5 | 0.8 | 0.3 | 0.25 | 0.9005519 | 0.0025916 | 0.0518327 | 0.7900100 | 0.9820180 | 154.5557517 | 0.0000000 |
| 5 | 0.8 | 0.3 | 1.00 | 0.5884416 | 0.0087466 | 0.1749316 | 0.2246004 | 0.8841159 | 10.1115560 | 0.0000000 |
| 5 | 0.8 | 0.3 | 4.00 | 0.1503721 | 0.0071233 | 0.1424664 | 0.0009990 | 0.4257493 | -49.0821650 | 0.0000000 |
| 5 | 0.8 | 1.0 | 0.25 | 0.7920380 | 0.0024577 | 0.0491536 | 0.7022727 | 0.8971029 | 118.8266235 | 0.0000000 |
| 5 | 0.8 | 1.0 | 1.00 | 0.5995904 | 0.0040889 | 0.0817777 | 0.4374875 | 0.7643856 | 24.3563692 | 0.0000000 |
| 5 | 0.8 | 1.0 | 4.00 | 0.1726299 | 0.0044209 | 0.0884190 | 0.0239011 | 0.3627373 | -74.0497284 | 0.0000000 |
| 5 | 1 | 0.1 | 0.25 | 0.8940559 | 0.0048582 | 0.0971640 | 0.6282717 | 1.0000000 | 81.1114992 | 0.0000000 |
| 5 | 1 | 0.1 | 1.00 | 0.4758841 | 0.0137388 | 0.2747767 | 0.0069181 | 0.9960040 | -1.7553080 | 0.0792066 |
| 5 | 1 | 0.1 | 4.00 | 0.1965135 | 0.0129869 | 0.2597375 | 0.0009990 | 0.9639361 | -23.3687092 | 0.0000000 |
| 5 | 1 | 0.3 | 0.25 | 0.8604745 | 0.0033270 | 0.0665396 | 0.7150849 | 0.9680320 | 108.3489209 | 0.0000000 |
| 5 | 1 | 0.3 | 1.00 | 0.4881194 | 0.0095941 | 0.1918820 | 0.1195804 | 0.8563686 | -1.2383254 | 0.2155954 |
| 5 | 1 | 0.3 | 4.00 | 0.1007742 | 0.0057848 | 0.1156956 | 0.0009990 | 0.4299451 | -69.0131201 | 0.0000000 |
| 5 | 1 | 1.0 | 0.25 | 0.7348826 | 0.0029112 | 0.0582234 | 0.6383616 | 0.8551698 | 80.6832908 | 0.0000000 |
| 5 | 1 | 1.0 | 1.00 | 0.4977522 | 0.0041455 | 0.0829107 | 0.3036963 | 0.6643606 | -0.5422101 | 0.5876738 |
| 5 | 1 | 1.0 | 4.00 | 0.1246628 | 0.0034988 | 0.0699754 | 0.0259241 | 0.2897103 | -107.2769047 | 0.0000000 |
| 5 | 1.25 | 0.1 | 0.25 | 0.8516084 | 0.0062154 | 0.1243090 | 0.5332667 | 0.9980020 | 56.5700727 | 0.0000000 |
| 5 | 1.25 | 0.1 | 1.00 | 0.4237737 | 0.0129779 | 0.2595572 | 0.0108891 | 0.9451049 | -5.8735634 | 0.0000000 |
| 5 | 1.25 | 0.1 | 4.00 | 0.1613786 | 0.0128102 | 0.2562041 | 0.0009990 | 0.9990260 | -26.4337184 | 0.0000000 |
| 5 | 1.25 | 0.3 | 0.25 | 0.8178372 | 0.0041144 | 0.0822890 | 0.6378871 | 0.9641359 | 77.2490361 | 0.0000000 |
| 5 | 1.25 | 0.3 | 1.00 | 0.4056618 | 0.0086395 | 0.1727904 | 0.1228521 | 0.7922328 | -10.9193729 | 0.0000000 |
| 5 | 1.25 | 0.3 | 4.00 | 0.0920504 | 0.0051176 | 0.1023529 | 0.0009990 | 0.3383866 | -79.7142890 | 0.0000000 |
| 5 | 1.25 | 1.0 | 0.25 | 0.6606144 | 0.0032006 | 0.0640114 | 0.5593906 | 0.8052448 | 50.1830249 | 0.0000000 |
| 5 | 1.25 | 1.0 | 1.00 | 0.4090260 | 0.0040690 | 0.0813807 | 0.2657343 | 0.5804446 | -22.3576409 | 0.0000000 |
| 5 | 1.25 | 1.0 | 4.00 | 0.0917183 | 0.0027186 | 0.0543716 | 0.0119630 | 0.2268981 | -150.1820280 | 0.0000000 |
5.3 Some fun
# Precompute some summarize grid
z <- out %>%
filter(b != 0) %>%
filter(scale == 1) %>%
group_by(b, rss, k) %>%
dplyr::select(on_toxic) %>%
mutate(
rss = as.numeric(as.character(rss)),
id = seq_along(on_toxic)
)
# Observed draws
a <- epred_draws(m, expand.grid("cat_pre_wt_log_scale" = c(-2, 2),
"beta" = c(-5, 5),
"var_trt" = c("low_var","high_var"),
"session_id" = "foo")) %>%
group_by(cat_pre_wt_log_scale, beta, draw) %>%
summarise(val = mean(val)) %>%
group_by(
cat_pre_wt_log_scale, draw
) %>%
summarise(
delta = diff(val)
) %>%
group_split(cat_pre_wt_log_scale, .keep = FALSE) %>%
lapply(
function(x){
x$delta
}
)
names(a) <- c("small", "large")
# Some simple wrapper for hypothesis testing later in the model scenarios
compute_delta <- function(hypothesis, res, ref_list){
temp <- hypothesis %>%
left_join(res, by = c("k", "rss","b"), relationship = "many-to-many") %>%
group_by(cat_size, id) %>%
summarise(
delta = diff(on_toxic)
)
out1 <- temp %>%
group_by(cat_size) %>%
do(as.data.frame(t(summarise_vec(.$delta)))) %>%
arrange(rev(cat_size))
out2 <- temp %>%
group_split(cat_pre_wt_log_scale, .keep = FALSE) %>%
suppressMessages() %>%
suppressWarnings() %>%
lapply(function(x){
x$delta
})
names(out2) <- c("large", "small")
r1 <- bayestestR::overlap(ref_list$small, out2$small)
r2 <- bayestestR::overlap(ref_list$large, out2$large)
return(
list(
"effects" = out1,
"overlap_coef" = sigfig(c("small" = r1, "large" = r2))
)
)
}5.4 Model scenarios
# No arrestment, no immigration
data.frame("k" = c(1, 1, 1, 1),
"rss" = c(1, 1, 1, 1),
"cat_size" = c("s", "s", "l", "l"),
"b" = c(-5, 5, -5, 5)) %>%
compute_delta(z, a) $effects
# A tibble: 2 × 6
# Groups: cat_size [2]
cat_size mean median sd lower upper
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 s -0.00205 0.00300 0.0891 -0.176 0.183
2 l -0.00205 0.00300 0.0891 -0.176 0.183
$overlap_coef
small large
"0.23" "0.13" # Size dependent arrestment, no immigration
data.frame("k" = c(0.25, 0.25, 4, 4),
"rss" = c(1,1, 1, 1),
"cat_size" = c("s", "s", "l", "l"),
"b" = c(-5, 5, -5, 5)) %>%
compute_delta(z, a)$effects
# A tibble: 2 × 6
# Groups: cat_size [2]
cat_size mean median sd lower upper
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 s 0.138 0.138 0.0639 0.0189 0.261
2 l -0.0983 -0.107 0.0803 -0.217 0.0783
$overlap_coef
small large
"0.44" "0.28" # No arrestment, size dependent immigration
data.frame("k" = c(1, 1, 1, 1),
"rss" = c(0.8,0.8, 1.25, 1.25),
"cat_size" = c("s", "s", "l", "l"),
"b" = c(-5, 5, -5, 5)) %>%
compute_delta(z, a)$effects
# A tibble: 2 × 6
# Groups: cat_size [2]
cat_size mean median sd lower upper
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 s 0.0190 0.0210 0.0880 -0.166 0.192
2 l -0.0110 -0.00849 0.0865 -0.185 0.155
$overlap_coef
small large
"0.26" "0.15" # Size dependent arrestment, size dependent immigration
data.frame("k" = c(0.25, 0.25, 4, 4),
"rss" = c(0.8,0.8, 1.25, 1.25),
"cat_size" = c("s", "s", "l", "l"),
"b" = c(-5, 5, -5, 5)) %>%
compute_delta(z, a)$effects
# A tibble: 2 × 6
# Groups: cat_size [2]
cat_size mean median sd lower upper
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 s 0.129 0.125 0.0568 0.0339 0.247
2 l -0.0753 -0.0804 0.0651 -0.188 0.0850
$overlap_coef
small large
"0.38" "0.18" # Size dependent arrestment in clustered treatment, no immigration
data.frame("k" = c(1, 0.25, 1, 4),
"rss" = c(1,1, 1, 1),
"cat_size" = c("s", "s", "l", "l"),
"b" = c(-5, 5, -5, 5)) %>%
compute_delta(z, a)$effects
# A tibble: 2 × 6
# Groups: cat_size [2]
cat_size mean median sd lower upper
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 s 0.235 0.235 0.0693 0.104 0.374
2 l -0.375 -0.380 0.0790 -0.500 -0.194
$overlap_coef
small large
"0.62" "0.58" # Size dependent arrestment in clustered treatment, size dependent immigration
data.frame("k" = c(1, 0.25, 1, 4),
"rss" = c(0.8,0.8, 1.25, 1.25),
"cat_size" = c("s", "s", "l", "l"),
"b" = c(-5, 5, -5, 5)) %>%
compute_delta(z, a)$effects
# A tibble: 2 × 6
# Groups: cat_size [2]
cat_size mean median sd lower upper
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 s 0.211 0.208 0.0592 0.110 0.331
2 l -0.328 -0.332 0.0655 -0.445 -0.166
$overlap_coef
small large
"0.53" "0.68" 5.5 Plots
5.5.1 Simplified plot
out <- read_csv("simulation/move_rules_sim.csv")
out <- out %>%
mutate(
scale = round(scale / 1000 * 12, 2),
rss = as.factor(round(exp(rss), 2))
)
out %>%
filter(b != 0) %>%
filter(scale == 1) %>%
group_by(b, rss, k) %>%
dplyr::select(on_toxic) %>%
mutate(
id = seq_along(on_toxic)
) %>%
arrange(b, rss, k, id) %>%
group_by(rss, k, id) %>%
summarise(
delta = diff(on_toxic)
) %>%
ggplot(aes(x = as.factor(k), y = delta, color = factor(rss))) +
geom_hline(aes(yintercept = 0), color = "grey",
size = 1, linetype = "dashed") +
geom_point(position = position_jitterdodge(jitter.width = 0.2,
jitter.height = 0,
dodge.width = 0.5),
alpha = 0.2,
size = 2) +
geom_point(stat = "summary",
position = position_dodge(width = 0.5),
shape = "—",
size = 4,
color = "black",
aes(group = factor(rss))) +
theme_bw(base_size = 15) +
theme(legend.position = "top") +
guides(colour = guide_legend(
override.aes = list(alpha = 1, size = 4),
title.position = "top"
)) +
labs(x = "Less toxic diet arresetment strength (ratio)",
y = expression(atop(Change~"in"~proprtion~time~on, paste("more toxic diet",~(beta[5]-beta[-5])))),
color = "Less toxic diet relative immigration strength (odds ratio)") +
scale_color_brewer(type = "seq") -> g;g
5.5.2 Full result
out <- read_csv("simulation/move_rules_sim.csv")
out <- out %>%
mutate(
scale = round(scale / 1000 * 12, 2),
rss = as.factor(round(exp(rss), 2))
)
k_lab <- unique(out$k)
names(k_lab) = sprintf("k = %s", k_lab)
scale_lab <- unique(out$scale)
names(scale_lab) = sprintf("scale = %s", scale_lab)
rss_lab <- unique(out$rss)
names(rss_lab) = sprintf("rss = %s", rss_lab)
out %>%
ggplot(aes(x = as.factor(k), y = on_toxic, color = factor(b))) +
geom_point(position = position_jitterdodge(jitter.width = 0.2,
jitter.height = 0,
dodge.width = 0.5),
alpha = 0.2,
size = 1) +
geom_point(stat = "summary",
position = position_dodge(width = 0.5),
shape = "—",
size = 3,
color = "black",
aes(group = factor(b))) +
theme_bw(base_size = 15) +
guides(colour = guide_legend(
override.aes = list(alpha = 1, size = 3),
)) +
facet_grid(scale ~ rss, labeller = labeller(
rss = reverse_names(rss_lab), scale = reverse_names(scale_lab)
)) +
guides(colour = guide_legend(override.aes = list(alpha = 1))) +
labs(x = "Less toxic diet arresetment strength (ratio)",
y = "Proportion time on toxic diet",
color = expression(beta)) +
scale_color_brewer(type = "qual")->g;g
6 Session Info
sessionInfo()R version 4.3.1 (2023-06-16 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)
Matrix products: default
locale:
[1] LC_COLLATE=English_United States.utf8
[2] LC_CTYPE=English_United States.utf8
[3] LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.utf8
time zone: America/New_York
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] spat1f_0.0.1 herbivar_0.2.1 forcats_0.5.1 stringr_1.5.0
[5] dplyr_1.1.2 purrr_1.0.1 readr_2.1.2 tidyr_1.3.0
[9] tibble_3.2.1 tidyverse_1.3.1 extraDistr_1.9.1 imager_0.42.13
[13] magrittr_2.0.3 doSNOW_1.0.20 snow_0.4-4 iterators_1.0.14
[17] foreach_1.5.2 glmmTMB_1.1.7 tidybayes_3.0.2 performance_0.10.4
[21] sjPlot_2.8.15 ggpubr_0.6.0 ggplot2_3.5.1 piecewiseSEM_2.3.0
[25] survival_3.3-1
loaded via a namespace (and not attached):
[1] fs_1.5.2 matrixStats_1.2.0 bitops_1.0-7
[4] sf_1.0-8 devtools_2.4.5 lubridate_1.8.0
[7] httr_1.4.3 RColorBrewer_1.1-3 insight_0.19.3
[10] doParallel_1.0.17 numDeriv_2016.8-1.1 profvis_0.3.7
[13] tools_4.3.1 backports_1.4.1 DT_0.23
[16] sjlabelled_1.2.0 utf8_1.2.2 R6_2.5.1
[19] mgcv_1.8-40 ggdist_3.1.1 urlchecker_1.0.1
[22] withr_2.5.0 sp_1.5-0 readbitmap_0.1.5
[25] gridExtra_2.3 Brobdingnag_1.2-9 prettyunits_1.1.1
[28] cli_3.6.1 shinyjs_2.1.0 sandwich_3.0-2
[31] labeling_0.4.2 sass_0.4.1 mvtnorm_1.1-3
[34] spatstat.data_3.0-0 brms_2.19.0 proxy_0.4-27
[37] StanHeaders_2.26.13 DHARMa_0.4.6 MuMIn_1.47.5
[40] sessioninfo_1.2.2 readxl_1.4.0 rstudioapi_0.13
[43] circular_0.5-0 visNetwork_2.1.0 generics_0.1.2
[46] gtools_3.9.2.2 crosstalk_1.2.0 vroom_1.5.7
[49] car_3.1-2 distributional_0.3.0 inline_0.3.19
[52] loo_2.5.1 Matrix_1.5-4.1 fansi_1.0.3
[55] abind_1.4-5 lifecycle_1.0.3 multcomp_1.4-25
[58] yaml_2.3.5 carData_3.0-5 grid_4.3.1
[61] qgam_1.3.4 promises_1.2.0.1 crayon_1.5.1
[64] miniUI_0.1.1.1 lattice_0.21-8 haven_2.5.0
[67] cowplot_1.1.1 pillar_1.9.0 knitr_1.43
[70] boot_1.3-28.1 estimability_1.3 shinystan_2.6.0
[73] codetools_0.2-19 imagerExtra_2.0.0 glue_1.6.2
[76] RcppArmadillo_0.11.2.0.0 V8_4.2.0 remotes_2.4.2
[79] vctrs_0.6.3 png_0.1-7 Rdpack_2.3.1
[82] cellranger_1.1.0 gtable_0.3.0 assertthat_0.2.1
[85] cachem_1.0.6 xfun_0.39 rbibutils_2.2.8
[88] mime_0.12 coda_0.19-4 shinythemes_1.2.0
[91] units_0.8-0 DiagrammeR_1.0.9 ellipsis_0.3.2
[94] fitdistrplus_1.1-11 TH.data_1.1-1 gap_1.2.3-6
[97] nlme_3.1-162 CircStats_0.2-6 xts_0.13.1
[100] usethis_2.1.6 bit64_4.0.5 threejs_0.3.3
[103] rstan_2.26.13 rprojroot_2.0.3 tensorA_0.36.2
[106] bslib_0.3.1 TMB_1.9.4 KernSmooth_2.23-20
[109] colorspace_2.0-3 DBI_1.1.3 tidyselect_1.2.0
[112] processx_3.6.1 emmeans_1.7.5 curl_4.3.2
[115] bit_4.0.4 compiler_4.3.1 amt_0.2.1.0
[118] rvest_1.0.2 arrayhelpers_1.1-0 see_0.7.1
[121] xml2_1.3.3 desc_1.4.1 colourpicker_1.1.1
[124] bayestestR_0.13.1 posterior_1.2.2 dygraphs_1.1.1.6
[127] checkmate_2.1.0 scales_1.3.0 classInt_0.4-7
[130] callr_3.7.0 tiff_0.1-11 digest_0.6.29
[133] fftwtools_0.9-11 spatstat.utils_3.0-1 minqa_1.2.4
[136] rmarkdown_2.24 base64enc_0.1-3 htmltools_0.5.2
[139] moveHMM_1.9 pkgconfig_2.0.3 jpeg_0.1-9
[142] lme4_1.1-29 highr_0.9 dbplyr_2.2.0
[145] fastmap_1.1.0 rlang_1.1.1 htmlwidgets_1.5.4
[148] shiny_1.7.1 farver_2.1.0 jquerylib_0.1.4
[151] svUnit_1.0.6 zoo_1.8-12 jsonlite_1.8.8
[154] bayesplot_1.10.0 geosphere_1.5-18 munsell_0.5.0
[157] Rcpp_1.0.8.3 stringi_1.7.12 MASS_7.3-60
[160] plyr_1.8.7 pkgbuild_1.3.1 parallel_4.3.1
[163] sjmisc_2.8.9 deldir_1.0-6 ggeffects_1.2.3
[166] splines_4.3.1 hms_1.1.1 sjstats_0.18.2
[169] ps_1.7.1 igraph_1.3.2 spatstat.geom_3.0-3
[172] markdown_1.1 ggsignif_0.6.3 reshape2_1.4.4
[175] stats4_4.3.1 pkgload_1.3.2 rstantools_2.2.0
[178] reprex_2.0.1 evaluate_0.15 RcppParallel_5.1.5
[181] modelr_0.1.8 bmp_0.3 nloptr_2.0.3
[184] tzdb_0.3.0 httpuv_1.6.5 RgoogleMaps_1.4.5.3
[187] polyclip_1.10-0 broom_1.0.5 gap.datasets_0.0.5
[190] xtable_1.8-4 e1071_1.7-11 rstatix_0.7.2
[193] later_1.3.0 class_7.3-22 memoise_2.0.1
[196] ggmap_3.0.2 bridgesampling_1.1-2